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Abstract 

There  are  many  situations  in  which  GPS  is  either  unable  to  provide  the  desired 
level  of  accuracy  or  is  unavailable.  Use  of  a  pseudolite-based  reference  system  for  nav¬ 
igation  can  be  a  means  for  positioning  during  these  times.  While  there  are  advantages 
in  using  a  pseudolite-based  reference  system,  there  are  still  implementation  issues  and 
deficiencies  that  must  be  addressed.  In  many  cases  a  pseudolite  system  with  ground 
based  transmitters  has  difficulty  in  determining  the  height  of  the  receiver  accurately. 
This  is  due  to  the  poor  vertical  observability  inherent  in  the  geometry  of  the  system. 
A  common  approach  in  naval  applications  for  solving  the  problem  of  poor  vertical 
observability  is  to  use  a  height  constraint,  which  is  well  known  when  travelling  on  a 
surface  of  water.  For  a  ground-based  vehicle,  knowledge  of  the  surface  topography 
can  be  obtained,  but  cannot  be  readily  used  in  the  same  manner  as  marine  cases, 
since  the  height  is  often  a  varying  function  of  position. 

This  research  investigates  and  develops  five  methods  of  incorporating  the  known 
surface  topography  in  a  non-linear  batch  least  squares  estimation  algorithm  using 
carrier-phase  measurements  from  pseudolites.  The  floating  point  carrier-phase  ambi¬ 
guities  are  estimated  in  this  process.  Real  and  simulated  data  sets  are  used  to  evaluate 
the  performance  of  the  five  algorithms.  In  simulation,  all  methods  performed  equiv¬ 
alently  well  on  a  flat  surface.  When  simulating  a  hill,  constraining  the  solution  to 
lie  in  a  plane  tangent  to  the  surface  topography  appeared  to  aid  the  solution  with 
the  best  knowledge  of  the  terrain.  Use  of  a  pseudo- measurement,  a  commonly  used 
approach,  did  not  provide  the  best  results,  and  indicates  the  inadequacy  of  using 
this  method  for  pseudolite-based  systems.  Results  using  data  from  a  real  system  on 
a  ground-based  vehicle  demonstrated  sub-decimeter  level  positioning  accuracy  in  all 
three  dimensions. 
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Methods  For  Aiding  F[eight  Determination 
In  Pseudolite-Based  Reference  Systems 
Using  Batch  Least-Squares  Estimation 

I.  Introduction 

1.1  Motivation 

The  Global  Positioning  System  (GPS)  has  become  an  unparalleled  aid  in  navi¬ 
gation.  GPS  was  created  as  a  military  program  with  limited  civilian  use.  Since 
becoming  operational  in  1993,  civilian  use  has  grown  rapidly  throught  the  world.  Not 
just  used  for  positioning,  GPS  gives  accurate  timing  capabilities  which  has  made  it  a 
cornerstone  of  countless  civil  systems. 

With  the  positioning  now  commonplace  with  GPS,  it  has  become  a  fight  to 
squeeze  as  much  capability  from  GPS  as  possible.  Research  in  both  the  civilian  and 
military  sectors  are  attempting  to  meet  the  addiction  of  accurate  positioning.  Ac¬ 
curacy  is  addicting.  Once  accurate  positions  to  meter  level  has  been  attained,  it 
is  unthinkable  to  stop  there.  Countless  applications  are  now  being  considered  that 
require  increasingly  higher  precision.  Many  systems  requiring  high  accuracy  often  re¬ 
quire  robustness  as  well.  Availability  of  GPS,  at  first  blush,  seems  all  encompassing. 
It  is  all  over  the  world.  But  looking  closer,  GPS  is  much  more  fragile  than  initially 
thought.  Accuracy  of  solutions  with  an  un-obstructed  sky  are  still  highly  dependent 
on  the  satellite  geometry.  With  any  occlusion  of  the  sky  by  obstacles,  there  can  be 
significant  loss  of  needed  accuracy.  Research  has  found  many  answers  to  these  prob¬ 
lems.  Differential  GPS  (DGPS),  carrier-phase  ambiguity  resolution,  and  augmented 
GPS  have  all  boosted  the  accuracy  and  robustness  of  GPS  based  navigation.  Many 
forms  of  DGPS  now  exist  to  aid  the  accuracy  of  GPS  navigation.  Integrated  GPS 
systems  primarily  combine  GPS  with  an  inertial  navigation  system  (INS)  and  create 


1 


a  more  robust  system.  Despite  these  improvements  new  demands  require  even  more 
of  GPS. 

A  flight  navigation  system  demands  some  of  the  highest  requirements  in  navi¬ 
gation.  The  United  States  Air  Force  (USAF)  is  interested  in  applications  that  need 
accuracy  within  tens  of  centimeters  or  less.  This  accuracy  must  be  maintained  with 
a  very  high  level  of  robustness.  A  loss  in  accuracy  for  even  a  short  time  could  be 
catastrophic  in  a  flight  system.  This  brings  the  problem  of  GPS  jamming  into  view. 
While  GPS  could  provide  the  desired  accuracy,  it  cannot  provide  the  robustness  to 
jamming  that  is  required. 

Development  of  systems  operating  in  GPS  jamming  environments  poses  another 
problem.  How  are  these  systems  tested?  A  truth  reference  must  have  an  accuracy  at 
least  an  order  of  magnitude  better  than  the  system  under  test.  In  the  past  systems 
were  tested  while  jamming  a  single  GPS  frequency,  allowing  the  truth  reference  system 
to  use  the  second  available  frequency  for  its  positioning.  In  reality  the  enemy  will  jam 
both  GPS  frequencies.  This  requires  that  new  systems  be  tested  in  this  environment. 
It  also  requires  that  the  flight  test  reference  system  also  function  under  dual  frequency 
jamming. 

The  746th  Test  Squadron’s  Central  Inertial  Guidance  Test  Facility  (CIGTF), 
Holloman  Air  Force  Base  (AFB),  NM  tests  and  evaluates  new  flight  navigation  sys¬ 
tems.  The  current  flight  reference  system  at  the  Holloman  AFB  incorporates  an 
Embedded  GPS/INS  (EGI),  a  second  GPS  receiver,  a  second  INS,  and  a  transpon¬ 
der/interrogator.  This  system  is  called  the  CIGTF  Reference  System  (CRS)  [14].  In 
the  past,  the  test  facility  was  able  to  test  the  robustness  of  GPS  based  flight  refer¬ 
ence  systems  in  the  presence  of  single  frequency  jamming.  This  testing  was  good  for 
research  purposes,  but  again  does  not  test  the  situation  that  will  most  likely  occur 
in  reality.  Currently  the  746th  does  not  have  the  capability  to  maintain  a  highly 
accurate  (near  centimeter  level)  truth  reference  under  dual-frequency  GPS  jamming 
conditions. 
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A  new  flight  test  reference  system  is  needed  that  can  meet  the  demands  of 
the  future.  This  reference  system  must  not  rely  on  GPS,  so  that  systems  designed 
to  operate  in  dual  frequency  jamming  can  be  tested  adequately.  Furthermore,  the 
accuracy  of  the  system  must  be  at  least  on  the  order  of  centimeters.  A  new  technology 
holds  the  answer  to  creating  this  reference  system. 

Pseudo-Satellites,  also  called  pseudolites,  are  transmitters  that  generate  a  GPS- 
like  signal  [10].  The  new  reference  system  will  be  a  pseudolite  based  reference  system. 
This  system  will  use  only  pseudolites  to  create  a  GPS-like  environment.  Being  able  to 
tailor  the  characteristics  of  the  pseudolites  will  allow  the  reference  system  to  operate 
on  a  frequency  other  than  GPS.  Signals  used  by  the  truth  reference  system  will  now  be 
unaffected  by  the  jamming  used  to  simulate  hostile  conditions  for  the  system  under 
test.  Using  pseudolites  will  solve  the  problem,  but  will  naturally  entail  new  issues 
that  research  must  overcome. 

Previous  research  has  simulated  the  environment  of  a  pseudolite  based  reference 
system  [8].  The  next  phase  of  research  is  creation  of  a  real  system.  Implementation  of 
a  proof-of- concept  is  underway  using  road  tests.  While  road  tests  are  only  a  stepping 
stone  to  creating  a  flight  reference  system,  which  inspired  this  research,  there  is  a 
large  amount  of  ground  based  positioning  that  requires  high  accuracy  and  the  ability 
to  function  in  a  hostile  environment.  Thus  solving  the  problems  of  ground  based 
pseudolite  positioning  is  extremely  important. 

1.2  Problem  Definition 

There  are  inherent  problems  that  arise  when  using  ground  based  pseudolites. 
The  greatest  of  these  lies  in  the  geometry  deficiencies.  With  all  the  transmitters  near 
the  ground  there  is  low  observability  in  the  vertical  direction.  This  must  be  overcome 
to  gain  three-dimensional  accuracies  on  the  order  of  centimeters. 

Suppose  the  position  of  a  static  receiver  is  being  estimated  using  pseudoranges. 
The  known  location  of  the  receiver  will  be  at  the  origin  of  some  ENU  reference  co- 


3 


Table  1.1:  Locations  of  pseudolites  in  ENU  frame 
with  the  receiver  at  the  origin.  The  slant  angle  gives 
the  angle  of  the  pseudolite  from  the  horizontal  plane 
as  seen  from  the  receiver. 


PRN  # 

East  (m) 

North  (m) 

Up  (m) 

Slant  Angle  (deg) 

1 

-1000 

1000 

100 

4.04 

2 

1000 

-1000 

10 

.405 

3 

-1000 

-1000 

150 

6.05 

4 

1000 

1000 

5 

.203 

ordinate  frame.  The  locations  of  4  pseudolites  in  the  same  ENU  frame  are  shown  in 
Table  1.1.  The  slant  angles  show  that  the  pseudolites  are  all  low  on  the  horizon.  This 
is  a  simple  case  that  demonstrates  poor  vertical  observability.  Using  least-squares 
estimation  with  single  differenced  pseudoranges  the  position  of  the  receiver  can  be 
found  with  some  uncertainty. 

To  simulate  this  example  the  true  ranges  were  calculated  and  used  as  the  pseu¬ 
dorange  measurements.  Measurement  error  was  introduced  as  a  simple  zero-mean 
white,  Gaussian  noise  on  the  measurements.  A  normally  distributed  random  noise 
with  strength  of  1  m2  was  added  to  each  measurement  (this  value  was  arbitrarily 
chosen  for  illustration  purposes).  These  measurements  were  then  used  to  estimate 
the  receiver  position  using  a  simple  ILS  method.  Estimation  is  done  in  3-dimensions. 
This  was  repeated  1000  times  to  show  the  distribution  of  the  solutions  obtained.  Fig¬ 
ures  1.1  and  1.2  show  the  results.  Figure  1.1  displays  the  1000  estimated  receiver 
positions.  The  scatter  plot  of  solution  points  has  been  projected  onto  the  3  axes  for 
easy  interpretation  of  the  errors  in  each  direction.  The  vertical  errors  can  be  hundreds 
of  meters  off,  despite  a  small  amount  of  measurement  noise.  Figure  1.2  highlights  on 
the  horizontal  errors.  The  east  errors  range  approximately  from  -10  to  10  meters  and 
the  north  errors  from  approximately  -4  to  4  meters. 

Having  poor  observability  in  the  vertical  dimension  increases  the  horizontal 
errors  due  to  the  non-linear  nature  of  this  problem.  Improving  the  vertical  error  will 
therefore  improve  the  horizontal  errors  as  well.  If  knowledge  of  the  height  is  known 
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Figure  1.1:  Receiver  position  estimation  error  in  3  dimensions. 
1000  runs  are  shown  with  projections  onto  each  axis  (shown  in 
red).  No  aiding  has  been  done. 
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Figure  1.2:  Horizontal  position  estimation  errors  with  no  aid¬ 
ing. 
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Figure  1.3:  Horizontal  position  estimation  errors.  Aiding  has 
been  done. 

from  a  source  other  than  the  given  measurements  this  can  be  used  to  aid  the  solution. 
The  simplest  method  is  to  assume  the  height  is  known.  Constraining  the  height  will 
allow  the  solution  to  be  estimated  in  the  horizontal  plane  only. 

This  was  simulated  and  the  results  are  shown  in  Figure  1.3.  Comparison  to  the 
horizontal  errors  in  Figure  1.2  will  show  a  reduction  of  errors  when  height  constraining 
is  used.  The  east  and  north  errors  have  both  been  reduced  to  ranges  of  -2  to  2  meters. 

The  previous  example  illustrates  what  simply  constraining  the  solution  to  a 
known  height  can  accomplish.  This  technique  is  very  simple  and  is  used  in  practice. 
Many  applications  allow  the  user  to  assume  a  receiver  height. 

In  a  ground  based  test  the  receiver  height  is  not  as  intrinsically  known  is  in 
some  cases.  The  receiver  is  also  in  motion.  This  does  not  mean  that  information  is 
not  available.  If  the  ground  topography  were  known  prior  to  the  test,  then  knowledge 
of  where  the  receiver  is  in  the  east  and  north  directions  will  indicate  the  receiver 
height.  Unfortunately,  each  direction  is  being  estimated  simultaneously.  Methods  of 
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incorporating  this  height  information  into  the  estimation  of  the  position  of  a  moving 
receiver  must  be  developed. 

This  research  focuses  on  evaluating  different  methods  of  incorporating  the  knowl¬ 
edge  of  the  height  into  the  solution.  The  height  information  is  obtained  from  assuming 
that  the  receiver  is  travelling  upon  an  unchanging  surface.  This  surface  can  be  sur¬ 
veyed  prior  to  testing,  and  then  used  to  aid  the  solution  algorithm. 

1.3  Proposed  Solution 

The  algorithms  that  will  be  developed  are  all  based  upon  using  least  squares 
estimation  (LSE).  Actually,  a  nonlinear,  iterative,  weighted,  batch  least  squares  esti¬ 
mation  method  will  be  implemented.  This  will  address  the  issues  in  dealing  with  a 
position  solution  using  carrier-phase  measurements  and  the  inherent  ambiguities. 

A  surface  grid  will  be  created  by  using  GPS  to  map  points  along  the  desired 
test  trajectory.  It  is  important  to  note  that  the  GPS  surveying  of  these  points  must 
also  have  centimeter  level  accuracy.  The  creation  of  the  grid  will  not  be  affected  by 
the  GPS  jamming  in  the  test,  since  the  grid  will  be  formed  prior  to  the  test.  In  the 
road  test  conducted,  no  system  was  being  tested  under  any  GPS  jamming,  therefore 
DGPS  measurements  were  collected  during  the  pseudolite  testing. 

Several  methods  for  integrating  the  information  from  the  grid  into  the  solution 
are  proposed  in  this  research: 

•  Method  1:  Only  the  horizontal  components  of  the  position  are  estimated.  The 
vertical  height  is  assumed  to  be  a  known  value.  Once  the  horizontal  position 
is  known,  it  is  used  to  find  the  corresponding  vertical  height  on  the  grid.  This 
new  height  is  used  in  the  next  iteration  as  the  height  constraint.  The  solution 
is  then  iterated  upon  in  this  fashion.  Thus  the  iterative  solution  is  constrained 
to  the  surface  grid. 

•  Method  2:  Instead  of  forcing  the  solution  to  be  in  the  horizontal  plane,  all  three 
dimensions  are  esimated.  The  low  observability  in  the  vertical  direction  that 
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would  normally  cause  the  solution  to  be  unobtainable  is  corrected  by  creating  a 
pseudo-measurement  using  the  GPS  grid.  This  measurement  is  created  much  the 
same  way  as  in  Method  1  (by  finding  the  corresponding  height  using  the  current 
horizontal  position  of  the  solution).  This  height  is  then  treated  as  an  additional 
measurement  and  is  incorporated  into  the  algorithm.  Since  a  weighted  LSE 
method  is  being  used,  the  pseudo-measurement  can  be  weighted  very  high  to 
keep  the  solution  very  close  to  the  actual  grid. 

Method  3:  Instead  of  assuming  that  the  2  horizontal  directions  are  the  best 
choices  to  constrain  the  system,  it  is  possible  to  find  the  2  best  axes.  By  using 
singular  value  decomposition  (SVD),  the  optimal  dimension  reduction  can  be 
made.  Depending  on  the  geometry  of  the  pseudolites  to  the  receiver,  the  two 
principal  axis  may  be  slightly  different  than  the  horizontal  directions  that  are 
normally  used.  A  state  transformation  is  used  so  that  the  solution  is  found  in 
the  optimal  plane  based  on  the  geometry.  There  will  still  be  problems  with 
the  vertical  direction.  The  SVD  approach  will  only  allow  the  solution  space  to 
be  optimized,  removing  the  direction  with  the  least  information;  however,  the 
solution  is  still  only  in  2-dimensions.  The  information  from  the  grid  will  be 
incorporated  by  constraining  the  solution  to  lie  on  the  known  surface.  Unlike 
Method  1  however,  the  solution  will  have  been  constrained  to  a  plane  that 
optimally  incorporates  the  measurements.  This  method  will  incorporate  the 
SVD  approach  at  each  epoch  when  forming  the  batch  process. 

Method  4:  This  Method  will  also  incorporate  the  SVD  approach.  It  will  com¬ 
prise  the  same  method  as  in  method  3  but  will  be  used  on  the  full  batch  process. 
This  will  allow  all  the  information  known  in  the  batch  process  to  aid  in  finding 
the  optimal  solution  space. 

Method  5:  Another  method  similar  to  the  SVD  approach  will  rely  upon  the 
geometry  of  the  GPS  grid  instead  of  the  pseudolite/receiver  geometry.  Instead  of 
finding  a  solution  on  a  plane  based  upon  the  relative  geometry  of  the  pseudolites 
and  receiver,  the  solution  will  find  the  plane  that  is  tangent  to  the  surface  of 


the  pre-defined  surface  at  the  current  position.  In  an  environment  that  has  a 
more  dynamic  surface  than  a  flat  plane,  this  method  will  allow  the  solution  to 
be  solved  in  the  directions  that  correspond  more  closely  to  the  known  height 
surface. 

1.4  Goals 

The  goal  of  this  research  is  to  implement  and  analyze  the  above  mentioned 
methods  for  aiding  a  pseudo lite-based  reference  system.  Since  real  data  is  available, 
it  will  be  used  in  the  analysis.  Simulation  data  will  also  be  generated  so  that  different 
surfaces  and  noise  sources  can  be  tested.  Each  method  will  be  created  and  tuned  to 
give  its  best  performance  with  real  data.  The  simulations  will  then  be  created  to  test 
each  algorithm’s  robustness  and  compare  them  against  each  other. 

1 . 5  Scope 

This  research  will  further  develop  the  algorithms  used  in  finding  the  navigation 
solution  using  a  pseudolite-based  reference  system.  The  real  data  that  has  been 
collected  is  from  a  pseudolite  test  facility.  The  data  has  been  collected  in  a  moving 
ground  vehicle.  Truth  trajectories  were  collected  at  the  same  time  as  the  pseudolite 
data.  The  algorithms  and  simulations  will  be  created  in  the  Mat  lab®  7.0  environment. 

Along  with  developing  the  algorithms  mentioned  above,  methods  to  deal  with 
real  data  issues  are  created.  Methods  to  find  carrier-phase  cycle  slips  are  addressed, 
since  there  are  different  problems  associated  with  pseudolites  than  with  normal  GPS 
receivers.  Furthermore,  methods  for  dealing  with  tropospheric  error  are  also  evalu¬ 
ated.  As  will  be  discussed  later,  tropospheric  error  must  be  handled  differently  than 
in  most  differential  GPS  cases. 

The  pseudolites  used  for  this  research  are  manufactured  by  Locata  Inc.,  and 
provide  a  time  synchronization  protocol  that  solves  the  problem  of  pseudolite  clock 
errors. 
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Figure  1.4:  The  Inverted  Pseudolite  Architecture. 

1.6  Assumptions 

The  following  assumptions  are  made  in  this  thesis: 

1.  Real-time  solutions  are  not  required,  because  the  research  is  directed  toward 
the  use  as  the  truth  trajectory  for  a  flight  test  reference  system.  All  data  is  able 
to  be  post  processed. 

2.  The  pseudolites  are  synchronized  using  the  manufacturer’s  TimeLoc  system. 
This  eliminates  the  time  offset  that  is  usually  inherent  with  pseudolites. 

3.  No  jamming  analysis  is  required,  because  the  pseudolites  operate  at  a  different 
frequency  than  GPS. 

1.7  Related  Research 

1.7.1  Inverted  Pseudolite  System.  In  April  1995,  the  746th  Test  Squadron 
implemented  a  novel  pseudolite  based  reference  system.  This  system  used  an  inverted 
architecture  shown  in  Figure  1.4.  In  the  test,  a  van  was  used  to  carry  a  pseudolite 
transmitter.  Another  pseudolite  was  placed  on  the  ground  as  a  reference  signal. 
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With  receivers  placed  at  known  locations  on  the  ground,  it  is  possible  to  implement  a 
reference  system.  This  research  was  a  proof-of-concept  for  creating  a  psendolite  based 
reference  system.  It  was  shown  to  have  accuracy  at  the  decimeter  level,  with  promise 
of  even  better  accuracy  with  further  development.  A  drawback  is  that  the  mobile  unit 
does  not  have  access  to  the  position  solution  unless  it  is  received  as  another  signal 
from  the  ground  station.  Further  detail  can  be  obtained  in  [22].  A  key  finding  of  the 
research  was  that  an  accurate  solution  is  highly  dependent  on  accurate  surveying  of 
the  reference  receivers  and  reference  pseudolite. 

1.7.2  Analysis  of  Pseudolite  Augmentation  for  GPS  Airborne  Applications. 

In  2002  the  Satellite  Navigation  and  Positioning  Group  (SNAP)  presented  research 
investigating  the  augmentation  of  pseudolites  with  GPS  for  airborne  applications  [15]. 
Implementation  issues  such  as  antenna  locations  and  placement  of  the  pseudolites  was 
considered.  Placement  of  a  second  antenna  on  the  bottom  of  the  aircraft  was  needed 
for  reception  of  ground-based  transmitters.  It  was  found  that  the  magnitude  of  the 
antenna  offset  (between  the  GPS  and  pseudolite  antennas)  will  determine  the  accuracy 
in  the  attitude  measurements.  Pseudolites  were  found  to  significantly  improve  the 
transmitter  geometry  and  the  signal  availability  of  a  airborne  navigation  system.  It 
was  also  found  that  the  addition  of  more  than  3  psedolites  into  the  system  gained 
only  marginal  improvement.  The  research  does  not  investigate  the  deficiencies  when 
landing,  although  in  this  research  both  pseudolites  and  GPS  satellites  are  used. 

1.7.3  Simulation  of  a  Pseudolite  Reference  System.  Development  and  simu¬ 
lation  of  a  pseudolite-based  reference  system  is  given  in  [8,9].  The  master’s  thesis  uses 
the  flight  path  of  a  C-12  as  simulation  truth  data.  Double  differenced  measurements 
are  used  as  well  as  optimal  smoothing  with  Kalman  filters.  The  tropospheric  error 
is  first  reduced  using  a  modified  troposphere  model,  and  then  the  residual  error  is 
reduced  through  the  differencing  techniques.  The  final  tropospheric  residual  error  is 
still  significant  enough  to  be  of  concern.  This  research  compares  two  different  methods 
to  remove  the  remaining  error.  One  is  a  weighted  covariance  method,  while  the  other 
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solution  is  to  specifically  model  the  error  by  adding  states  to  the  Kalman  filter.  This 
research  showed  that  single  frequency  operation  is  sufficient  for  centimeter  level  accu¬ 
racy.  This  means  that  widelane  techniques  are  not  needed.  Another  finding  confirmed 
that  a  change  in  relative  geometry  is  required  to  resolve  the  carrier-phase  ambiguities. 
Finally,  this  research  showed  the  dominant  errors  in  a  simulated  pseudolite  system  are 
from  residual  tropospheric  delay  and  misalignment  of  the  pseudolites.  This  research 
lacks  insights  into  real  problems  that  may  occur  with  actual  measurements. 

1.7.4  Naval  Height  Constraining .  ffeight  constraining  is  very  applicable  in 

naval  position  estimation  as  the  height  of  a  body  of  water  is  often  well  known.  [27]. 

1.7-4-1  Improvement  of  Ambiguity  Resolution  with  Constraints.  In 
this  research,  [27],  decreasing  the  time  it  takes  to  resolve  ambiguities  on  the  fly  on¬ 
board  a  ship  with  two  roving  receivers  was  investigated.  Several  constraints  were 
used  in  an  attempt  to  aid  the  ambiguity  resolution  process.  A  height  constraint  was 
used  as  well  as  a  heading  and  pitch  constraint,  a  fixed  baseline  between  the  rover 
pair,  and  an  ambiguity  constraint.  A  pseudo- measurement  was  used  to  incorporate 
the  known  height.  The  height  constraint  was  found  to  be  the  most  beneficial  of  the 
constraints  used  in  reducing  the  resolution  time  by  nearly  half.  It  was  also  found  that 
the  reliability  of  the  system  was  increased  from  80  to  95  percent. 

1.7. 4-2  Navigation  in  Constricted  Waterways.  Constricted  waterways 
demand  precise  positions  for  safe  navigation.  This  is  difficult  given  that  many  times 
the  line-of-sight  GPS  signals  are  blocked  by  obstructions.  Two  research  techniques  are 
described  below  to  improve  navigation  in  this  environment.  Both  researchers  include 
a  height  constraint  incorporated  as  a  pseudo-measurement. 

Research  by  [17]  involves  augmenting  GPS  with  pseudolites  to  give  more  ac¬ 
curacy  and  robustness  in  cases  when  GPS  is  masked  by  obstructions.  It  was  found 
that  using  even  a  single  pseudolite  the  HDOP  is  decreased  dramatically  as  well  as  the 
worst-case  horizontal  error. 
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Augmenting  GPS  with  GLONASS  and  a  height  constraint  is  a  second  method  of 
improving  the  accuracy  and  robustness.  Using  the  second  constellation  of  positioning 
satellites  results  in  an  increased  number  of  visible  transmitters  at  any  given  time.  The 
addition  of  a  height  constraint  makes  a  further  improvement  to  the  reliability,  since 
the  height  is  usually  known  for  marine  applications  [23]. 

1.7.5  LocataLite  Test  Examples.  LocataLites  have  been  used  in  several 
experimental  trials  as  pilot  studies  of  the  system. 

1.7. 5.1  BlueScope  Steel.  In  this  pilot  study,  [2],  the  Locata  system  is 
shown  to  give  centimeter  level  relative  positioning  precision  in  a  severe  multipath  en¬ 
vironment.  The  BlueScope  steelworks  at  Port  Kernbla  in  Australia  uses  a  large  crane 
to  move  slabs  of  steel.  The  position  of  the  crane  must  be  tracked  to  assure  accurate 
placement  and  inventory  inside  a  large  warehouse.  Due  to  the  metallic  structure  the 
environment  is  highly  conducive  to  extreme  multipath.  Nevertheless,  the  Locata  sys¬ 
tem  was  able  to  give  sub-centimeter  accurate  positioning  of  the  crane  compared  to  a 
total  station.  Multipath  was  mitigated  by  careful  placement  of  LocataLites  and  Lo¬ 
cata  receivers.  Future  testing  will  also  give  the  steelworks  the  ability  to  track  vehicles 
and  steel  slabs  moving  from  outdoors  to  indoor  environments. 

1.7. 5. 2  Structural  Deformation.  In  another  study,  LocataLites  were 
used  to  measure  the  structural  deformation  of  a  suspension  bridge  [3].  The  test 
occurred  at  the  suspension  footbridge  at  Parsley  Bay  at  Sydney  Harbor,  Australia. 
Accurate  measuring  of  bridge  deformation  can  be  difficult  to  attain  with  GPS  systems, 
since  their  accuracy  can  fluctuate  significantly  depending  on  satellite  geometry  and 
availability.  A  Locata  receiver  was  placed  on  the  bridge  with  four  LocataLites  placed 
at  lower  altitudes.  The  system  was  initialized  by  manually  surveying  without  the  use 
of  GPS.  It  was  shown  that  LocataLites  gave  sub-centimeter  level  relative  positioning 
accuracy  without  the  use  of  GPS  signals. 
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1.8  Thesis  Overview 


Chapter  2  presents  all  the  pertinent  background  theory  needed  in  this  research. 
The  least  squares  estimation  methods,  GPS  fundamentals,  and  pseudolite  character¬ 
istics  will  be  discussed  in-depth.  Chapter  3  details  the  methods  developed  in  this 
research  for  aiding  the  navigation  solution,  detecting  cycle-slips,  and  creating  sim¬ 
ulations.  Chapter  4  analyzes  the  results  of  the  simulations  and  the  processing  of 
the  real  data,  as  well  as  compare  the  different  methods  with  one  another.  Chapter 
5  summarizes  the  results  and  also  provides  recommendations  for  future  research  in 
boosting  the  accuracy  of  the  different  methods  for  aiding  the  position  solution  of  a 
pseudolite-based  reference  system. 
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II.  Background 


2. 1  Overview 

First,  the  GPS  signal  structure  and  measurements,  which  are  similar  to  paseu- 
dolite  measurements,  will  be  discussed.  The  next  sections  will  introduce  least 
squares  estimation  methods.  These  will  form  the  basis  of  the  methods  developed  in 
later  solution  algorithms.  After  discussing  least  squares  estimation,  pseudolites  are 
introduced,  and  they  are  compared  and  contrasted  with  GPS.  The  errors  and  issues 
that  corresponds  to  pseudolites  are  then  detailed.  Finally,  a  short  introduction  to  the 
specific  pseudolites,  from  Locata  Corporation,  is  made. 

2.2  GPS  Signals 

In-depth  presentation  of  GPS  signals  can  be  found  in  [16,19].  GPS  signals  oper¬ 
ate  fundamentally  on  two  frequencies.  These  two  center  frequencies  are  1575.42  Mhz 
and  1227.6  Mhz,  termed  LI  and  L2  respectively.  The  LI  channel  transmits  three  mes¬ 
sages,  a  coarse-acquisition  (C/A)  code,  a  precise  (P)  code,  and  a  navigation  message. 
L2  only  contains  the  P  code  and  the  navigation  message.  The  signals  are  transmitted 
on  the  carrier  frequency  using  a  bi-phase  shift  key  (BPSK)  modulation.  The  GPS 
satellites  have  very  expensive  clocks  onboard  that  guarantee  the  signals  are  trans¬ 
mitted  at  synchronized  and  precise  times.  With  knowledge  of  the  satellite’s  location, 
given  through  ephcmeris  data  in  the  navigation  message,  and  knowledge  of  the  signal 
structure,  a  receiver  can  track  a  satellite’s  signal  and  knowing  the  transmission  time 
can  create  a  pseudorange  measurement.  This  measurement  is  not  the  true  range  to 
the  satellite  due  to  the  effects  of  clock  errors.  Another  measurement  that  is  available 
to  the  receiver  is  a  carrier-phase  measurement.  A  receiver  tracks  the  carrier  signals 
phase  component  and  can  determine  a  very  accurate  measure  of  distance  to  the  satel¬ 
lite  as  long  as  the  number  of  complete  signal  cycles  that  occurred  before  reception  is 
discovered.  The  phase  can  be  measured  very  accurately,  which  can  be  converted  to 
a  distance  using  the  wavelength  of  the  signal.  Thus,  with  a  wavelength  of  19.03  cm 
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using  the  LI  signal  and  24.42  cm  using  the  L2  signal,  a  carrier-phase  measurement 
can  be  accurate  to  within  centimeters  by  measuring  phase  within  l/20th  of  a  cycle. 

Pseudolite  signals  are  designed  very  similar  to  GPS.  Similar  methods  for  attain¬ 
ing  pseudorange  and  carrier-phase  measurements  are  used.  A  third  measurement  that 
will  be  exploited  is  the  signal  power  measurement.  The  signal-to- noise  ratio  (SNR) 
will  be  measured  at  the  receiver  to  give  indications  of  poor  or  bad  measurements. 

2.2.1  Pseudorange  Measurements.  Pseudorange,  commonly  called  code, 
measurements  can  be  expressed  as: 

p  =  r  +  c(5tu  -  5tsv)  +T  +  /  +  mp  +  vp  (2.1) 


where 

p  =  pseudorange  measurement  (meters) 
r  =  true  range  between  receiver  and  transmitter  (meters) 
c  =  speed  of  light  (meters  /  second) 

Stu  =  receiver  clock  error  (seconds) 

5tsv  =  transmitter  clock  error  (seconds) 

T  =  errors  due  to  troposphere  (meters) 

/  =  errors  due  to  ionosphere  (meters) 

Trip  =  errors  due  to  pseudorange  multipath  (meters) 
vp  =  pseudorange  errors  due  to  receiver  noise  (meters) 

The  time  errors  that  are  associated  with  this  measurement  come  directly  from 
the  errors  in  the  clocks  in  the  receiver  and  the  transmitter.  For  GPS,  the  transmitter 
clock  error  is  very  small,  since  atomic  clocks  are  being  used.  The  receiver  clock  is 
usually  not  an  atomic  clock,  and  thus  contains  large  errors.  This  receiver  clock  error 
must  either  be  estimated  or  eliminated  through  differencing. 

Nearly  the  same  measurement  model  is  used  with  pseudolites,  except  for  the 
ionospheric  error  term.  Since  pseudolites  are  generally  ground  based,  they  do  not 
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travel  through  the  ionosphere,  and  thus  do  not  contain  ionospheric  errors.  It  is  im¬ 
portant  to  note  that  with  pseudolites  the  transmitter  clocks  will  not  be  atomic  based, 
and  therefore  the  transmitter  clock  error  may  be  much  larger  than  with  GPS.  As  will 
be  seen  later,  the  pseudolite  used  in  this  research  handles  the  problem  of  transmitter 
clock  error  in  a  unique  way. 

2.2.2  Carrier-Phase  Measurements.  Carrier-phase  measurements  are  simi¬ 
lar  to  code  measurements,  but  they  can  be  much  more  accurate  since  the  phase  of  a 
signal  can  be  measured  with  very  high  accuracy.  This  measurement  can  be  expressed 
as 


0  =  A  1(r  +  c(5tu  -  Stpi)  +  T  —  I  +  +  v^)  +  N  (2.2) 

where 

0  =  Carrier-phase  measurement  (cycles) 

A  =  carrier  wavelength  (meters  /  cycle) 
r  =  true  range  between  receiver  and  transmitter  (meters) 
c  =  speed  of  light  (meters/second) 

5tu  =  receiver  clock  error  (seconds) 

Stpi  =  transmiter  clock  error  (seconds) 

T  =  errors  due  to  troposphere  (meters) 

/  =  errors  due  to  ionosphere  (meters) 

rricj,  =  errors  due  to  carrier-phase  multipath  (meters) 

V'f,  =  carrier-phase  errors  due  to  receiver  noise  (meters) 

N  =  carrier-phase  integer  ambiguity  (cycles) 

Note  that  the  sign  of  the  ionospheric  error  is  reversed  from  the  code  case.  This 
occurs  because  the  phase  of  the  carrier  is  advanced  in  the  ionosphere  instead  of  de¬ 
layed,  as  occurs  in  the  pseudorange.  This  creates  the  problem  of  code-carrier  di¬ 
vergence  in  the  presence  of  large  ionospheric  error.  Again  the  ionospheric  delay  is 
neglected  when  using  the  above  measurements  with  pseudolites. 
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The  carrier-phase  ambiguity  term  is  unknown  and  must  be  resolved  before  an 
accurate  solution  is  obtained.  This  process  is  called  ’’ambiguity  resolution,”  and  it  is 
described  in  the  next  section. 

2.2.3  Carrier-Phase  Ambiguity  Resolution.  When  using  carrier-phase  mea¬ 
surements,  the  integer  ambiguity  term  must  be  resolved  to  achieve  accurate  position¬ 
ing.  The  integer  ambiguity  represents  the  unknown  number  of  cycles  that  are  present 
before  integration  of  the  Doppler  measurement  begins.  The  ambiguities  are  constant 
for  the  measurement  history  as  long  as  the  measurements  are  continuous  and  free 
of  cycle  slips.  Solving  for  the  ambiguities  usually  requires  enough  angular  motion 
between  the  reference  system  and  the  mobile  receiver  for  the  correct  ambiguities  to 
become  apparent.  In  GPS,  this  can  take  tens  of  minutes  since  the  angular  motion 
between  a  ground  based  mobile  receiver  and  the  satellite  constellation  is  slow.  In  a 
pseudolite  system,  the  angular  motion  can  be  much  more  active  and  the  time  required 
for  resolution  of  the  ambiguity  may  be  much  shorter. 

The  ambiguity  is  in  reality  an  integer,  since  it  represents  a  measure  of  whole 
cycles.  It  is  difficult  to  force  an  algorithm  to  estimate  the  ambiguity  state  as  an 
integer.  This  is  doubly  difficult  since  any  errors  that  are  bias  like  in  the  system  can 
be  indistinguishable  from  being  part  of  the  ambiguity.  This  gives  rise  to  a  method  of 
fixing  the  ambiguities  to  integers  after  they  have  been  estimated  as  floating  values. 

There  are  many  ways  that  the  integer  ambiguity  is  found,  many  rely  on  the  use 
of  Kalman  filter  estimates  of  the  ambiguities  and  the  covariances.  Examples  of  this 
are  the  Least  squares  AMbiguity  Decorrelation  Algorithm  (LAMBDA)  and  the  Fast 
Ambiguity  Search  Filter  (FASF).  Both  methods  are  described  in  [20]. 

In  this  research  the  ambiguities  will  remain  floating.  This  means  that  integer 
ambiguity  resolution  will  not  be  done.  In  the  current  stage  of  development,  the  Locata 
system  outputs  phase  measurements  that  have  non-integer  ambiguities. 
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Figure  2.1:  Single  Difference  Between  Receiver  A  and  Pseu- 
dolites  1  and  2 

2.2.4  Single  Differencing.  Single  differencing  is  used  to  eliminate  or  reduce 
some  errors  in  measurements.  When  linear  combinations  of  measurements  between 
receivers,  transmitters,  or  times  are  taken  the  common  errors  can  be  removed.  Single 
differencing  between  two  transmitters  (pseudolites)  will  use  the  notation  (V). 

Figure  2.1  shows  how  single  differencing  between  two  pseudolites  would  appear. 
Single  differencing  with  code  measurements  as  in  Figure  2.1  is  defined  as 

Vpa  =  Pa-Pa  (2-3) 

where  p\  is  the  code  measurement  between  receiver  A  and  pseudolite  1,  and  p\  is  the 
code  measurement  between  receiver  A  and  pseudolite  2. 

Single  differencing  will  eliminate  the  receiver  clock  error  and  reduces  the  tropo¬ 
spheric  errors.  This  can  be  seen  by  combining  Equations  (2.1)  and  (2.3)  to  yield 
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(2.4) 


Combining  like  terms  gives 


(2.5) 


The  receiver  clock  errors  5th  and  8ii  are  the  same  for  receiver  A  so  they  will  cancel. 

UA  UA  J 

The  remaining  differences  will  be  represented  as  (V).  Equation  (2.5)  can  now  be 
rewritten  as 


Vp\2  =  Vri2  +  cV*#  +  VT]2  +  VmJ,2  +  Vu 


,12 

Pa 


(2.6) 


In  the  GPS  case,  the  differencing  will  reduce  the  tropospheric  error,  but  will  increase 
the  multipath  and  measurement  noise  by  a  factor  of  y/2. 

Single  differencing  with  carrier-phase  measurements  is  very  similar  to  the  above 
code  differencing  and  yields 


V<^2  =  A_1(Vr^2  +  cVhfJ2  +  VTf  +  Vm“  +  VuJ2  )  +  VA^2  (2.7) 


The  receiver  clock  error  has  again  dropped  out,  and  the  differenced  integer  ambiguity 
term  has  appeared.  Examples  of  double  differencing  are  shown  in  [8,16,19,20].  This 
technique  is  used  to  eliminate  both  receiver  and  pseudolite  clock  errors.  It  also  allows 
further  reduction  of  the  atmospheric  errors,  but  with  a  penalty  of  further  increased 
multipath  and  measurement  noise. 
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2.3  Least  Squares  Estimation 

Least  squares  estimation  is  a  very  useful  and  common  tool  in  estimation.  Ex¬ 
amples  and  similar  derivations  can  be  found  in  [7,16].  Given  measurements  of  some 
system,  it  is  possible  to  find  the  characteristics  of  that  system.  The  characteristics  of 
the  system  are  termed  the  system  states.  They  must  be  related  to  the  measurements 
in  some  known  way.  Suppose  there  is  a  linear  relationship: 

z  =  Hx  +  v,  (2.8) 

where  z  is  the  m-dimensional  measurement  vector,  H  is  the  m  x  n  measurement 
matrix,  v  is  a  vector  of  m  unknown  errors,  and  x  is  the  n-dimensional  unknown  state 
vector.  If  the  number  of  measurements,  m,  is  greater  than  the  number  of  states,  n, 
the  system  is  over-determined.  This  is  the  common  case  in  this  research.  If  m  where 
to  be  less  than  n,  the  system  would  be  under-determined,  and  no  solution  would  be 
found.  For  a  solution  to  be  determined,  H  must  also  have  a  rank  of  n.  If  this  is  not 
the  case,  then  there  is  a  problem  with  unobservability,  which  will  be  discussed  later. 
The  goal  is  to  find  the  best  estimate  of  the  states,  x,  given  the  measurements. 

The  LSE  method  attempts  to  minimize  the  difference  in  a  least-squares  sense 
between  the  measurements  and  an  estimate  of  the  measurements  based  upon  estimates 
of  the  unknown  states,  x  and  the  known  matrix  H.  The  least  squares  scaler  error  is 
written  as 

J  =  (z  -  Hx)t  (z  -  Hx)  (2.9) 

The  estimate  of  the  vector  x  is  found  by  minimizing  the  scalar  error  with  respect  to 
the  state  vector  x.  This  is  done  by  setting  the  derivative  of  J  with  respect  to  (wrt)  x 
to  zero,  and  solving  for  x.  Once  this  is  done,  the  solution  for  x  comes  forth  as 

x  =  (HtH)-1  Htz  (2.10) 
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If  H  is  not  of  full  rank  then  the  inverse  of  H 1  H  is  singular.  Problems  also  occur 
if  (HrH)  1  is  nearly  singular,  which  is  the  case  with  poor  observability.  In  the 
observable  case,  the  least  squares  estimate  of  the  states  is  obtained. 


2.3.1  Weighted  LSE.  When  trying  to  estimate  the  states  of  a  system  it  is 
important  to  use  all  the  information  possible.  When  considering  measurements  it  is 
important  to  incorporate  the  knowledge  of  the  measurement  errors.  Given  Equation 
(2.8),  the  measurement  error  vector  can  be  assumed  to  be  zero-mean  white  Gaussian 
noise  with  covariance  matrix  R,  where 


R  =  £{vvt} 


(2.11) 


The  measremtent  noise  covariance  matrix  will  effectively  represent  the  errors  that  are 
expected  in  the  measurements  and  also  any  correlation  of  errors  between  the  different 
measurements. 

The  noise  covariance  matrix  is  applied  in  least  squares  estimation  and  forms 
the  method  commonly  called  weighted  least  squares  (WLS).  A  slight  change  occurs 
in  the  LSE  algorithm  and  now  appears  as 


(2.12) 


As  will  be  seen  later,  different  measurements  with  different  noise  variances  will  get 
weighted  differently  in  the  WLS  structure.  This  will  allow  more  accurate  measure¬ 
ments  to  have  more  influence  on  the  estimate. 

2-4  Non-Linear  Least  Squares  Estimaion 

Previously,  the  relationship  between  the  states  and  the  measurements  was  as¬ 
sumed  to  be  linear.  In  many  situations  this  is  not  a  valid  assumption.  Instead  of 
having  a  direct  linear  relationship,  the  measurements  will  be  a  nonlinear  representa¬ 
tion  of  the  system  states.  This  is  represented  below. 
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z  —  h  (x)  +  v 


(2.13) 


Since  h  (x)  is  potentially  nonlinear,  the  previous  method  of  least  squares  esti¬ 
mation  cannot  be  used.  The  iterative  least  squares  (ILS)  will  be  used.  This  method 
starts  with  an  initial  estimate  of  the  states  xG  and  uses  the  measurements  to  calculate 
a  corrected  state  estimate  xc.  The  corrected  state  estimate  then  becomes  the  current 
state  estimate  x,  which  replaces  the  initial  estimate.  This  process  is  repeated  until 
there  is  convergence.  Convergence  of  the  solution  is  tested  before  the  estimate  is 
validated.  Convergence  of  the  solution  can  depend  on  several  factors  such  as  accu¬ 
racy  of  the  initial  estimate  and  observability  of  the  states.  Details  of  this  process  are 
described  below. 

2-4-1  Iterative  Least  Squares.  The  iterative  estimation  method  involves 
correcting  the  estimates  of  the  states  at  each  iteration.  To  update  the  current  state 
estimate  and  create  the  corrected  state  estimate,  the  following  relationship  is  defined. 

xc  =  x  +  Ax  (2-14) 

where  Ax  is  defined  as  the  estimated  error  in  the  state  estimate. 

The  measurement  residuals  are  the  difference  between  the  actual  measurement 
vector  z,  and  the  modeled  measurements  given  the  current  state  estimate  x. 

Az  =  z  —  h  (x)  (2-15) 

The  measurement  residual  vector  Az  can  also  be  represented  using  the  state  correction 
Ax.  A  first-order  expansion  is  used  to  relate  the  measurement  residual  to  the  error 
in  the  state  estimate  as  follows: 


Az  =  HAx  +  v, 


(2.16) 
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where  Ax  is  the  error  in  the  state  estimate.  The  sensitivity  matrix  H  is  the  same 
as  above  and  embodies  the  relationship  of  the  states  and  measurements.  Now  the 
relationship  is  non-linear,  and  H  is  the  matrix  of  partial  derivatives  evaluated  from 
the  current  state  estimate. 


H 


dh 

<9x 


(2.17) 


The  Ah  row  of  H  is  made  up  of  the  partial  derivatives  of  the  zth  entry  in  the  h  vector 
wrt  the  states  x.  Therefore.  H  is  an  m  by  n  matrix,  with  m  remaining  the  number 
of  measurements  and  n  the  number  of  states. 

Since  second-order  terms  have  been  ignored,  the  previous  development  with  a 
linear  term  to  obtain  Equation  (2.10)  can  be  followed  to  similarly  obtain  the  following, 


Ax=(HtR_1H)  1  HrR  1  Az  (2.18) 

It  is  important  to  note  that  during  each  iteration  the  H  matrix  must  be  re-evaluated 
at  the  current  state  estimate.  Each  iteration  will  use  the  corrected  state  estimate 
from  the  previous  iteration  to  linearize  h  and  form  H. 

As  stated,  this  process  of  correcting  the  state  estimate  is  repeated  until  sufficient 
convergence  is  obtained.  Convergence  is  indicated  by  the  change  in  the  state  estimate 
(Ax)  becoming  sufficiently  small.  A  threshold  can  be  set  that  is  based  upon  knowledge 
of  the  system  and  measurements.  Using  a  sum  of  squares  of  Ax,  the  solution  is 
either  declared  valid  or  another  iteration  is  begun.  In  this  manner  the  non-linear 
state  estimate  is  obtained.  The  accuracy  of  the  above  linearization  is  dependent  on 
assuming  the  second-order  effects  of  the  states  in  the  system  are  minimal. 

2.5  Batch  Least-Squares 

The  above  methods  of  Least-Squares  Estimation  are  commonly  used  in  an  epoch 
by  epoch  process.  This  means  that  the  state  vector  and  measurements  for  each  time 
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epoch  are  used  independently  from  previous  and  subsequent  times.  The  states  at 
a  given  epoch  are  estimated  using  the  measurements  pertaining  to  that  time  before 
another  LSE  process  is  started  for  the  next  epoch.  This  is  common  to  the  techniques 
used  in  Kalman  Filter  applications. 

This  method  ignores  the  common  elements  between  epochs.  A  batch  process 
considers  all  the  states  and  measurements  together.  The  states  are  not  constrained  to 
be  related  to  any  single  epoch.  Everything  is  estimated  simultaneously.  This  allows, 
for  example,  the  ambiguity  states  that  are  common  to  many  epochs,  to  be  estimated 
properly. 

Estimating  all  the  states  using  all  the  measurements  can  involve  much  larger 
matrices  than  normally  encountered.  This  will  motivate  some  ad  hoc  techniques  to 
minimize  the  size  discussed  in  the  previous  chapter. 

The  batch  ILS  method  will  allow  other  errors  that  are  only  observable  through 
geometry  changes  as  well.  Errors  such  as  residual  tropospheric  delay  and  transmitter 
position  error  are  future  avenues  of  research  that  are  expected  to  see  benefit  using 
batch  ILS  estimation. 

2. 6  Pseudolites 

Pseudolites,  standing  for  pseudo-satellites,  are  usually  ground-based  transmit¬ 
ters  that  propagate  GPS-like  signals.  Pseudolites  have  been  used  to  augment  GPS 
for  many  years  [10,26].  One  reason  for  pseudolite  augmentation  of  GPS  is  to  boost 
the  availability  of  GPS  by  adding  ground  based  signals  that  are  copies  of  the  GPS 
signal.  This  allows  locations  with  poor  satellite  visibility  to  achieve  position  solutions. 
Pseudolites  are  much  cheaper  than  satellites  and  are  much  more  easily  customized. 
Different  algorithms  and  signal  types  can  be  used  to  mitigate  different  types  of  error. 
A  thorough  discussion  of  pseudolites  can  be  found  in  the  following  references  [10,26]. 

Many  of  the  errors  in  the  pseudolite  measurements  are  similar  to  the  errors  in 
the  GPS  measurements.  Tropospheric  error,  multipath,  receiver  noise,  etc.  are  all 
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similar.  Other  pseudolite  implementation  issues  are  synchronization  and  the  near/far 
problem.  Ionospheric  error  will  be  completely  neglected,  since  the  signals  will  not  be 
propagating  through  the  ionosphere.  These  error  sources  will  be  discussed  briefly,  and 
the  techniques  to  mitigate  them  will  be  discussed  after  the  measurements  themselves 
are  introduced. 

2.6.1  GPS  vs  Pseudolites.  The  core  difference  between  GPS  and  pseudolites 
lie  in  the  location  of  the  transmitters.  GPS  transmitters  are  located  on  satellites 
20,000  km  in  the  sky.  Pseudolites  can  be  placed  on  the  ground  in  almost  any  location. 
The  following  are  a  list  of  many  differences  between  GPS  and  pseudolites. 

•  Pseudolites  can  operate  on  frequencies  other  than  the  LI  and  L2  frequencies 
associated  with  GPS.  This  is  key  for  the  purposes  of  this  research,  since  the 
system  will  be  used  in  the  presence  of  jamming  on  these  two  frequencies. 

•  Pseudolite  signals  are  based  on  the  ground  and  therefore  do  not  pass  through 
the  ionosphere  as  GPS  signals  do.  Therefore  the  usual  ionospheric  errors  are 
not  present. 

•  Angular  motion  between  a  mobile  receiver  in  motion  and  the  transmitters  is 
greater  in  pseudolite  systems  then  with  GPS.  This  maybe  related  to  how  quickly 
carrier-phase  ambiguities  are  resolved. 

•  Multipath  in  a  psedolite  system  can  be  much  more  severe.  The  elevation  angles 
of  the  transmitters  can  be  much  lower  than  experienced  in  GPS.  With  much 
smaller  glance  angles  there  is  much  more  position  correlated  multipath. 

•  While  signal  strength  with  pseudolites  can  be  much  more  powerful  than  with 
GPS,  the  power  at  the  receiver  has  the  potential  to  be  highly  dynamic.  This 
can  cause  problems  in  receiver  design. 

•  Instead  of  orbital  erros  in  GPS,  pseudolite  systems  are  very  sensitive  to  accurate 
surveying  of  transmitter  position  errors.  The  physical  position  of  the  transmitter 
and  the  phase  center  of  the  antenna  must  be  found  accurately. 
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Figure  2.2:  Illustration  of  Near/Far  Problem 

2.6.2  Pseudolite  Errors.  The  following  sections  provide  further  discussion 
of  errors  that  occur  with  the  use  of  pseudolites  and  the  methods  used  to  mitigate 
these  errors. 

2.6.2. 1  Near-Far  Problem.  The  near/far  problem  is  an  issue  for  both 
pseudolite  augmented  GPS  systems  and  for  pseudolite  only  based  systems.  Maximiz¬ 
ing  the  transmitters’  useful  distance  requires  increase  in  the  transmission  levels.  The 
problem  occurs  as  a  receiver  approaches  a  pseudolite,  the  signal  from  that  pseudolite 
drowns  out  the  signals  from  the  other  pseudolites,  making  it  impossible  to  have  a  so¬ 
lution.  Figure  2.2  shows  the  donut  region  that  is  the  usable  area  around  a  pseudolite. 
A  higher  quality  pseudolite  will  have  a  larger  usable  region. 

The  near/far  problem  can  be  ignored  in  some  applications  if  the  usable  area 
is  acceptable,  but  this  is  not  typically  the  case.  The  standard  method  to  solve  this 
problem  is  to  use  a  pulsed  signal  scheme.  In  this  type  of  system,  each  pseudolite  will 
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transmit  on  a  predefined  duty  cycle.  This  allows  a  receiver  to  track  other  desired 
signals  while  in  close  proximity  to  a  given  pseudolite.  This  has  been  used  with  GPS 
augmentation  for  landing  systems  in  the  past.  The  problem  of  synchronization  now 
becomes  tied  to  the  near/far  problem,  because  each  individual  pseudolites  must  be 
synchronized  so  that  it  is  pulsing  at  the  correct  time  relative  to  the  other  pseudo¬ 
lites.  The  pulsing  scheme  is  commonly  referred  to  as  Time  Division  Multiple  Access 
(TDMA).  Other  forms  of  signal  separation  that  can  be  adopted  from  the  communi¬ 
cations  industry  are  Frequency  Division  Multiple  Access  (FDMA)  and  Code  Division 
Multiple  Access  (CDMA).  More  detailed  discussion  of  these  schemes  are  found  in  [10]. 

2. 6. 2. 2  Synchronization.  Synchronization  of  the  signals  is  vital  to  ref¬ 
erence  systems.  In  GPS,  the  satellites  use  extremely  accurate  atomic  clocks.  These 
clocks  are  very  expensive  and  not  feasible  for  pseudolites.  The  signals  must  be  syn¬ 
chronized  so  that  the  user  can  determine  what  time  it  was  transmitted.  This  is  the 
only  way  the  user  can  determine  his  pseudorange  to  the  pseudolite.  In  many  cases 
when  pseudolites  where  used  to  augment  GPS,  the  pseudolites  are  all  synchronized 
to  GPS  time  from  the  satellite  signals  themselves.  This  required  the  availability  of 
GPS  at  the  location  of  each  pseudolite.  GPS  will  be  unavailable  in  this  test,  and  the 
pseudolites  must  be  synchronized  using  a  different  method. 

To  overcome  this  problem  a  new  technology  has  been  developed  by  a  company 
called  Locata,  whereby  each  pseudolite  has  the  ability  to  precisely  synchronize  with 
other  pseudolites  in  the  system.  This  feature  and  others  will  be  discussed  later  in 
Section  2.7. 


2. 6. 2. 3  Tropospheric  Error.  A  large  source  of  errors  in  the  measure¬ 
ments  is  due  to  the  troposphere.  The  troposphere  is  the  portion  of  the  atmosphere 
that  is  neutral,  roughly  0-10km  in  altitude  [25].  In  the  proposed  system,  the  signals 
will  be  propagating  completely  in  the  troposphere.  The  delay  experienced  by  the  sig¬ 
nals  depend  on  the  refractivity  index  of  the  air  mass  the  signal  is  travelling  through. 
This  index  is  a  function  of  the  density  of  the  air,  which  is  in  turn  due  to  the  density 


of  the  dry  and  wet  components  of  the  air  [16].  Roughly  90%  of  the  delay  is  due  to 
the  dry  portion  of  the  density.  For  GPS  signals,  this  portion  can  be  predicted  to 
within  1%  @  zenith  (90  degree  elevation).  The  remaining  10%  due  to  the  water  vapor 
is  much  harder  to  calculate  and  is  only  predictable  to  10-20%  at  zenith  [19].  For 
pseudolite  applications,  the  errors  are  amplified.  A  commonly  used  model  with  GPS 
is  the  Hopfield  model.  It  is  based  on  predicting  the  refractivity  index  at  an  altitude 
based  on  the  refractivity  at  the  surface.  The  tropospheric  error  is  usually  calculated 
at  zenith  and  then  transformed  to  the  desired  elevation  using  a  mapping  function. 
These  methods  can  be  inadequate  for  pseudolite  applications,  because  the  mapping 
functions  become  inaccurate  at  lower  elevations,  which  occur  regularly  for  ground- 
based  pseudolites.  More  applicable  models  must  be  introduced  for  use  in  pseudolite 
applications. 

The  tropospheric  model  to  be  used  can  be  found  in  [8,9,11].  It  is  a  function 
of  temperature,  atmospheric  pressure,  relative  humidity,  elevation  angle,  and  range. 
The  atmospheric  parameters  are  measured  at  a  ground  based  reference  station.  The 
delay  for  the  mobile  receiver  is  defined  as 
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where  the  variables  are  defined  as 
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tapl,u  =  tropospheric  delay  for  mobile  receiver  (meters) 

Ru  =  slant  range  between  the  pseudolite  and  user  (meters) 
A hu  =  the  height  of  the  user  above  the  pseudolite  (meters) 
A rVtdry  =  differential  vertical  dry  delay  (meters) 

A TVjWet  =  differential  vertical  wet  delay  (meters) 


el, 


=  elevation  angle  in  radians 


A Hapl  =  difference  in  height  between  pseudolites  and  reference  receiver 
hs  =  height  of  reference  receiver 

Ps  =  surface  pressure  (millibars) 

Ts  =  surface  temperature  (Kelvins) 

Ns  =  surface  refractivity 

The  tropospheric  equation  is  valid  for  all  elevation  angles  except  zero.  A  modified 
equation  for  this  indeterminate  case  can  be  found  in  [11], 


After  the  majority  of  the  tropospheric  error  is  eliminated  through  use  of  the 
model  and  reduced  by  single  differencing,  the  remaining  error  is  termed  the  residual 
tropospheric  error.  Although  this  is  still  very  small,  it  cannot  be  ignored  when  seeking 
centimeter  level  accuracy.  This  residual  error  can  be  modelled  as  a  scale  factor  and 
reduced  or  incorporated  as  a  state  to  be  estimated.  These  methods  are  part  of  the 
research  in  this  thesis  and  are  discussed  later  in  this  paper.  Additional  details  can  be 
found  in  [8,9]. 


2. 6. 2-4  Multipath  Error.  Multipath  error  in  a  pseudolite-based  refer¬ 
ence  system  is  similar  to  using  GPS.  The  source  of  the  error  is  still  the  same.  Signals 
from  the  sending  antenna  travel  a  true  path  to  the  receiving  antenna,  but  they  can 
also  be  reflected  off  the  surrounding  landscape  to  the  receiving  antenna.  The  error 
is  caused  by  these  reflections  travelling  multiple  reception  paths.  In  pseudolites,  this 
error  is  often  amplified.  The  problem  is  increased  with  the  low  glance  angles  some¬ 
times  associated  with  pseudolite  locations.  A  more  complex  problem  lies  in  the  use 
of  differential  techniques  and  synchronizatiion  with  pseudolites. 
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In  GPS  applications  the  satellites  are  moving  and  a  stationary  receiver  will 
not  have  a  biased  multipath  since  the  source  of  the  signal  is  moving  throughout  the 
sky.  In  a  pseudolite  system  the  transmitters  do  not  move  relative  to  each  other. 
This  can  create  a  standing  multipath  dilemma  between  pseudolites.  In  a  differential 
pseudolite  system  the  standing  multipath  is  found  between  the  reference  receiver 
and  the  pseudolites.  In  the  test  designed  for  this  research,  the  pseudolites  use  a 
synchronization  scheme  requiring  the  pseudolites  to  receive  signals  from  each  other. 
This  leads  to  a  standing  multipath  problem. 

Multipath  error  can  be  difficult  to  eliminate.  In  GPS  applications,  special  an¬ 
tennas  can  be  used  to  reduce  most  multipath  errors.  This  is  easier  with  satellite 
applications,  since  most  usable  signals  are  all  above  15  degrees.  In  pseudolite  applica¬ 
tions,  signals  are  coming  from  ground  based  transmitters  and  force  different  receiver 
antennas. 


2. 6. 2. 5  Geometry/Observability.  In  GPS,  geometry  plays  a  large  role 
in  the  accuracy  of  the  solution.  Poor  geometry  can  easily  be  the  largest  source  of 
error.  Many  times,  waiting  for  better  satellite  geometry  is  the  only  solution.  This 
is  one  of  the  reasons  that  GPS  is  sometimes  augmented  with  pseudolites.  Having 
satellites  only  low  on  the  horizons  results  in  a  larger  error  in  the  vertical  direction. 
This  is  due  to  the  poor  observability  in  that  direction  by  the  measurements.  The 
same  problem  occurs  with  pseudolites. 

In  a  landing  scenario  with  ground  based  pseudolites,  at  the  moment  of  landing 
the  receiver  and  all  the  pseudolites  are  very  close  to  being  in  the  same  plane.  As 
they  come  closer  to  lying  in  the  same  plane,  the  ability  of  the  system  to  determine 
the  vertical  position  deteriorates.  This  could  be  solved  similarly  to  the  case  with 
GPS,  except  instead  of  waiting  for  a  satellite  to  move  into  a  better  position,  the 
pseudolites  can  be  placed  in  a  better  geometry  or  more  pseudolites  can  be  added 
to  the  existing  setup  to  create  a  better  geometry.  Both  approaches  have  the  same 
problem  though.  To  achieve  better  observability  in  the  vertical  direction,  some  of  the 
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pseudolites  must  be  placed  such  that  there  is  a  significant  elevation  angle  between  it 
and  the  receiver  (in  order  to  make  the  measurement  geometry  non-coplanar) .  This  can 
be  very  problematic  at  an  airfield,  where  tall  objects  and  landscapes  have  been  avoided 
for  safety  concerns.  Additionally,  placing  pseudolites  on  towers  maybe  undesirable, 
due  to  the  swaying  of  the  towers.  This  will  cause  a  very  significant  change  in  the 
location  of  the  transmitter  and  will  introduce  an  unacceptable  error,  when  centimeter 
accuracy  is  needed. 

Even  in  GPS  cases  geometry /observability  becomes  a  problem.  In  naval  navi¬ 
gation  there  can  be  cases  of  obstructed  views  of  satellites.  With  poor  visibility  the 
observability  of  the  vertical  dimension  can  be  poor.  In  many  naval  cases  the  height 
of  the  water  level,  thus  the  receiver  height,  is  well  known,  and  can  be  incorporated  as 
a  constraint  in  the  solution  [17,23,27].  Such  a  height  constraint  can  be  implemented 
by  using  a  pseudo-measurement  of  the  height. 

2. 6. 2. 6  Receiver  Noise.  Receivers  not  only  measure  the  desired  signals, 
but  also  introduce  unwanted  random  noise.  This  error  is  caused  by  many  sources, 
and  is  dependent  upon  receiver  design  and  location.  The  receiver  measurements  are 
affected  by  RF  radiation  in  the  band  of  interest;  the  antenna,  amplifiers,  cables,  and 
the  receiver;  multi-access  noise  (i.e.,  interference  from  other  GPS  signals  and  GPS- like 
broadcasts  from  system  augmentations);  and  signal  quantization  noise  [16].  Receiver 
noise  error  is  not  a  bias  and  can  be  mitigated  by  filtering. 

2. 7  LocataLites 

The  pseudolites  that  will  be  used  for  this  research  are  from  Locata  Corpora¬ 
tion  [1,4-6].  These  pseudolites,  termed  LocataLites,  are  a  novel  design  that  solves 
some  of  the  problems  discussed  in  the  previous  sections.  LocataLites  can  be  used  to 
set  up  a  LocataNet,  which  is  a  network  of  synchronized  LocataLites,  and  has  potential 
for  near-centimeter-level  accuracy  using  carrier-phase  positioning.  These  pseudolites 
solve  the  problem  of  synchronization  with  a  method  called  TimeLoc.  TimeLoc  syn- 
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chronizes  LocataLites,  and  can  be  described  in  a  series  of  steps  [4]. 


1.  LocataLite  A  transmits  a  C/A  code  and  carrier  signal  on  a  particular  PRN  code. 

2.  The  receiver  section  of  LocataLite  B  acquires,  tracks  and  measures  the  signal 
(C/A  code  and  carrier-phase  measurements)  generated  by  LocataLite  A. 

3.  LocataLite  B  generates  its  own  C/A  code  and  carrier  signal  on  a  different  PRN 
code  to  A. 

4.  LocataLite  B  calculates  the  difference  between  the  code  and  carrier  of  the  re¬ 
ceived  signal  and  its  own  locally  generated  signal.  Ignoring  propagation  errors, 
the  differences  between  the  two  signals  are  due  to  the  differences  in  the  clocks 
between  the  two  devices,  and  the  geometric  separation  between  them. 

5.  LocataLite  B  adjusts  its  local  oscillator  using  Direct  Digital  Synthesis  (DDS) 
technology  to  bring  the  code  and  carrier  differences  between  itself  and  LocataL¬ 
ite  A  to  zero.  The  code  and  carrier  differences  between  LocataLite  A  and  B 
are  continually  monitored  so  that  they  remain  zero.  In  other  words,  the  local 
oscillator  of  B  follows  precisely  that  of  A. 

6.  The  final  stage  is  to  correct  for  the  geometrical  offset  between  LocataLite  A  and 
B,  using  the  known  coordinates  of  the  LocataLites,  and  after  this  TimeLoc  is 
achieved. 

This  method  of  synchronization  can  be  replicated  with  many  LocataLites  creat¬ 
ing  a  network  of  synchronized  pseudolites  without  the  use  of  expensive  atomic  clocks. 
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III.  Algorithm  Development 


3. 1  Overview 

This  chapter  will  describe  the  algorithms  developed  for  aiding  a  pseudolite  ref¬ 
erence  system  with  height  determination.  Since  a  non-linear,  iterative,  batch 
least  squares  estimation  method  forms  the  basis  of  each  algorithm,  the  initial  setup 
for  each  algorithm  is  very  similar,  and  will  be  discussed  first.  Following  the  initial 
setup,  the  surface  grid  information  will  be  introduced.  An  important  aspect  of  height 
determination  methods  in  this  research  involves  the  pre-defined  surface  created  from 
GPS  surveys.  Finally,  the  algorithms  that  were  developed  will  be  introduced  next. 
Each  method  will  be  covered,  and  important  characteristics  of  each  will  be  illustrated. 

3.2  Initial  Algorithm  Setup 

3.2.1  Carrier  Phase  Ambiguities.  For  precise  positioning,  it  is  mandatory 
that  carrier-phase  measurements  be  used.  As  discussed,  this  means  that  the  related 
ambiguities  must  be  determined.  A  batch  least  squares  technique  is  used  in  this 
research  for  this  purpose. 

The  ambiguities  are  found  through  geometry  change.  Having  a  large  time  his¬ 
tory  of  measurements  without  any  motion  will  result  in  little/no  observability  into  the 
ambiguity  states.  Along  the  same  lines,  there  is  little  gain  in  having  a  time  history  of 
1000  epochs  over  that  of  10  epochs  if  there  is  the  same  amount  of  geometry  change. 
This  motivates  the  idea  of  pruning  the  epochs  that  will  be  used  in  the  ambiguity 
search  process.  Creating  a  batch  process  with  an  over-abundance  of  measurements 
that  will  not  gain  any  more  useful  information  only  impedes  the  solution  process.  The 
ILS  solution  contains  an  n  by  n  matrix  inversion,  where  n  is  the  number  of  states. 
Each  epoch  used  in  the  solution  adds  2  or  3  states,  depending  on  the  solution  type, 
to  the  overall  state  vector.  This  inversion  can  quickly  become  a  computational  prob¬ 
lem  for  a  processor.  Allowing  the  algorithms  to  only  select  out  a  predefined  number 
of  epochs  spanning  the  time  that  an  ambiguity  is  present,  captures  all  the  needed 
geometry  change  that  is  available  and  minimizes  the  number  of  states.  After  these 
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ambiguity  states  are  found,  an  epoch  by  epoch  ILS  position  solution  can  be  run  with 
the  known  ambiguities  no  longer  being  estimated.  This  final  solution  method  has  no 
need  to  be  a  batch  process  and  is  very  straightforward. 

One  major  problem  still  exist.  Using  the  batch  process,  each  ambiguity  state  is 
tied  to  the  measurements  at  each  epoch.  Any  time  there  is  a  cycle  slip,  the  ambiguity 
can  change.  From  that  time  forward  the  ambiguity  state  before  the  slip  is  no  longer 
associated  with  any  new  measurements.  This  requires  that  cycle  slips  be  identified.  If 
a  cycle  slip  goes  undetected,  the  ambiguity  estimate  that  is  attained  will  be  completely 
invalid,  and  will  cause  large  errors  in  the  position  solution. 

3. 2. 1.1  Epoch  Selection.  Assuming  cycle  slip  detection  has  been  done 
(which  is  discussed  in  the  following  section),  the  complete  set  of  ambiguities  and 
the  epochs  and  measurements  they  are  associated  with  are  known.  The  true  epochs 
that  will  be  used  in  the  batch  ILS  process  must  be  selected.  The  goal  of  the  batch 
process  is  to  find  the  floating  point  ambiguities,  so  it  must  select  the  epochs  that  span 
the  lifetime  of  all  the  ambiguities.  Each  transmitter,  and  its  related  carrier-phase 
measurements,  will  have  a  predetermined  number  of  ambiguity  states.  Selecting  a 
given  number  of  epochs  evenly  spaced  over  the  life  of  each  ambiguity  should  maximize 
the  observability  of  the  ambiguities  for  that  specific  transmitter.  This  approach  is 
taken  for  each  transmitter.  Since  the  number  of  ambiguities  and  their  individual  time 
periods  of  validity  do  not  necessarily  coincide  between  transmitters,  it  is  necessary  to 
include  all  the  needed  epochs  for  each  transmitter.  This  results  in  some  epochs  that 
are  duplicated  and  some  that  are  very  close  together.  Some  epochs  can  be  eliminated 
to  keep  the  size  of  the  batch  process  small.  Duplicated  times  can  first  be  eliminated, 
and  then  times  that  are  close  together  can  be  combined  by  merging  them  together  and 
using  an  epoch  between  the  two.  This  will  maintain  the  needed  epochs  for  estimation 
of  the  ambiguities,  while  creating  a  smaller  number  of  measurements  and  states  for  a 
manageable  batch  process. 
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Figure  3.1:  Example  of  epoch  selection  for  batch  processing. 

There  are  three  ambiguities  shown  that  span  different  times. 

Each  has  diamonds  indicating  preferred  epochs  spaced  evenly 
over  their  time  of  validity.  Along  the  x-axis  there  are  squares 
that  represent  times  chosen  for  use  in  the  batch  process. 

The  epoch  selection  and  pruning  process  is  illustrated  in  Figure  3.1.  Given 
three  ambiguities  that  span  the  time  history  with  differing  validity  regions,  epochs 
must  be  selected  that  represent  the  ambiguities  and  are  not  repeated  or  too  close 
together.  Each  ambiguity  has  initially  evenly  space  epochs  that  would  make  ideal 
choices  to  represent  that  ambiguity.  With  multiple  ambiguities  it  can  be  seen  that 
there  are  repeated  epochs  and  some  epochs  that  are  close  together.  Along  the  bottom 
of  the  figure  are  selected  epochs  that  are  projected  vertically  to  show  how  some  epoch 
possibilities  are  merged  to  form  a  single  epoch  choice.  In  this  case,  17  of  the  24 
unique  epoch  possibilities  are  selected  for  the  processing  of  the  batch  least  squares 
estimation. 


3.2.2  Cycle  Slip  Detection.  Cycle  slip  detection  is  necessary  when  processing 
carrier-phase  measurements.  In  some  receiver  designs  an  output  can  indicate  if  a 
cycle  slip  has  occurred.  Using  the  change  in  phase  measurements  between  epochs, 
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Figure  3.2:  Example  of  a  typical  SNR  plot  shown  with  bad 
measurements  indicated  as  being  below  an  SNR  bound.  In  this 
case  there  are  two  times  that  the  SNR  threshold  is  crossed. 

it  is  possible  to  indicate  a  cycle  slip  from  an  uncommonly  large  change  in  the  phase 
measurement  itself.  In  GPS  cases,  a  cycle  slip  happens  instantaneously,  making  it 
easier  to  identify  the  erroneous  event. 

In  its  current  configuration,  the  LocataLite  can  make  cycle  slip  detection  more 
difficult.  The  receiver  is  not  able  to  output  a  cycle  slip  indicator,  and  cycle  slips  do 
not  appear  instantaneously  as  in  GPS.  From  observation  of  test  data,  it  is  apparent 
cycle  slips  can  occur  over  a  gradual  time.  Although,  the  receiver  remains  in  frequency 
lock  with  the  measured  signal,  no  effort  is  made  to  maintain  the  phase  at  zero.  During 
a  period  of  high  multipath,  this  allows  the  phase  to  drift  slowly  and  cause  a  drifting 
cycle  slip  instead  of  an  instantaneous  error. 

In  an  attempt  to  locate  the  cycle  slips,  the  SNR  measurements  are  utilized  as 
an  indicator  of  a  weak  signal,  and  thus  potential  slips.  If  the  SNR  measurement  for  a 
given  transmitter  drops  below  a  pre-determined  threshold,  the  measurements  for  that 
transmitter  are  suspect.  More  specifically,  the  carrier-phase  measurement  is  assumed 
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to  potentially  have  a  cycle  slip.  This  form  of  detection  assumes  cycle  slips  are  only 
found  in  periods  of  low  SNR.  An  example  of  this  is  shown  in  Figure  3.2.  These  SNR 
measurements  are  made  at  the  receiver  and  correspond  with  a  single  transmitter.  This 
is  an  indication  of  the  received  power  of  the  signal  from  the  pseudolite.  In  the  figure, 
the  SNR  threshold  is  set  at  6  dB.  This  measurement  was  indicated  from  empirical 
evidence  from  the  Locata  Corporation.  There  are  two  measurements  that  fall  below 
this  bound.  For  the  carrier-phase  measurements  from  this  transmitter  there  would 
be  two  ambiguity  states.  One  would  be  associated  with  the  measurements  before  the 
possible  cycle  slip  and  the  other  would  be  associated  after  this  time. 

An  illustration  of  how  a  single  transmitter  can  have  several  ambiguity  states  is 
shown  in  Figure  3.3.  The  initial  ambiguity  related  with  this  transmitter  is  set  as  the 
first  ambiguity  state.  After  each  detected  cycle  slip,  there  must  be  a  new  ambiguity 
state  for  the  future  carrier  phase  measurements  from  this  transmitter.  This  can  be 
seen  as  the  ambiguity  state  changes  after  each  time  it  drops  to  0,  which  is  a  period 
of  invalid  measurements.  The  carrier-phase  measurements  cannot  be  used  at  these 
times.  After  each  possible  cycle  slip,  a  new  ambiguity  must  be  created.  This  explains 
why  an  ambiguity  state  is  never  reused.  It  is  important  to  remember  cycle  slips  are 
only  assumed  to  be  highly  possible  in  these  regions,  as  indicated  by  low  SNR  values. 

3.2.3  Initial  Trajectories.  While  carrier-phase  measurements  will  be  used  for 
precise  positioning,  the  code  measurements  will  not  be  thrown  away.  It  is  important 
in  a  batch  process  that  an  initial  trajectory  be  used  as  the  starting  point  for  the  ILS 
estimation.  In  a  normal  LSE  routine,  the  positions  are  found  epoch  by  epoch,  allowing 
the  solution  for  the  previous  epoch  to  be  used  as  the  initial  position  for  the  next  time. 
In  a  batch  process,  there  must  be  an  initial  guess  for  every  epoch  simultaneously. 
Furthermore,  this  initial  guess  must  be  different  for  at  least  some  epochs.  If  a  single 
position  were  used  for  all  epochs  as  the  initial  guess,  the  calculation  of  the  sensitivity 
matrix  H  would  produce  no  information  for  determining  the  ambiguities.  These 
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Figure  3.3:  Example  of  the  changing  ambiguity  states  for  a 
single  transmitter.  This  plot  shows  which  ambiguity  state  is 
associated  with  this  transmitter  at  each  epoch.  An  ambiguity 
state  of  0  indicates  the  SNR  from  this  transmitter  is  below  the 
threshold,  and  a  cycle  slip  has  likely  occurred. 

states  are  dependent  on  geometry  change,  starting  each  epoch  at  the  same  point 
causes  immediate  divergence. 

The  initial  trajectory  that  is  needed  for  formulation  of  the  batch  process  is  found 
using  the  pseudorange  measurements.  This  initial  trajectory  must  only  be  accurate 
enough  to  allow  the  batch  ILS  process  to  converge.  The  pseudorange  measurements 
are  generally  accurate  enough  for  this  purpose.  Since  single  differenced  measurements 
are  being  used,  there  is  no  state  associated  with  the  removed  receiver  clock  bias.  The 
initial  positions  are  calculated  using  the  same  ILS  process,  but  an  epoch  by  epoch 
approach  is  used  instead  of  a  batch  process. 

3.2.4  Frame  Rotations.  This  research  will  base  all  calculations  in  the  East 
North  Up  (ENU)  coordinate  frame.  Many  algorithms  only  focus  on  2-dimensional 
solutions.  Given  the  ground  based  nature  of  pseudolite  systems,  the  horizontal  plane 
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that  is  tangent  to  the  surface  of  the  earth  is  the  best  solution  plane.  Confining  the 
solution  to  this  plane  in  Earth  Centered  Earth  Fixed  (ECEF)  coordinates  would  be 
overly  difficult.  In  choosing  the  ENU  frame,  the  solution  can  be  confined  to  the 
horizontal  plane  and  elimination  of  the  vertical  direction  is  trivial.  Conversion  of  the 
pseudolite  locations  from  ECEF  to  ENU  coordinates  is  accomplished.  After  an  origin 
for  the  ENU  frame  is  chosen  all  the  locations  are  calculated  in  the  new  frame. 

A  further  coordinate  rotation  was  used  for  use  with  real  data.  In  the  real  world 
road  test  conducted  in  relation  to  this  research,  a  road  was  mapped  and  traversed 
using  GPS  receivers  and  a  Locata  receiver.  The  road  travelled  was  not  aligned  with 
East  or  North.  In  later  sections  the  creation  of  the  GPS  surveyed  grid  will  be  detailed, 
but  creating  the  grid  in  the  standard  ENU  frame  required  a  large  number  of  invalid 
grid  points  due  to  the  nature  of  the  road  geometry.  Instead  of  the  road  spanning  the 
length  of  a  rectangle,  which  is  desirable  for  a  larger  number  of  useful  grid  points,  the 
road  only  spanned  a  diagonal  section  of  a  more  square  grid.  To  facilitate  a  better  use 
of  computing  resources,  the  coordinate  frame  was  rotated  by  5  degrees.  This  aligned 
the  road  so  as  to  maximize  the  number  of  valid  points  in  the  grid. 

Figure  3.4  illustrates  this  concept  by  providing  the  exact  example  explained 
above.  Figure  (a)  shows  the  test  road,  viewed  from  above,  in  the  normal  ENU  frame. 
The  road  track  is  along  the  diagonal  of  where  the  grid  will  lie.  All  of  the  empty  space 
in  the  grid  will  be  invalid  points  that  cannot  be  used,  but  are  present  and  increase 
the  calculation  time.  Figure  (b)  shows  the  same  road,  also  viewed  from  above,  but 
this  time  the  ENU  frame  has  been  rotated  5  degrees,  so  the  road  now  fills  the  area 
where  the  grid  will  lie.  This  increases  the  number  of  useful  grid  points. 

3.2.5  Measurement  Formulation.  Measurements  from  each  pseudolite  are 
logged  in  the  receiver.  Pseudoranges,  carrier-phase,  and  SNR  measurements  are  avail¬ 
able.  As  previously  discussed,  single  differenced  measurements  will  be  used  in  order 
to  eliminate  the  receiver  clock  error.  A  base  transmitter  will  be  chosen  that  will  then 
be  differenced  with  the  remaining  transmitters.  This  creates  a  set  of  independent 
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(a)  (b) 

Figure  3.4:  (a)  Horizontal  view  of  a  test  road  in  ENU  frame.  Poor  use  of  grid. 

(b)  Horizontal  view  of  a  test  road  in  rotated  ENU  frame.  Road  now  spans  grid  more 
efficiently 

combinations  of  the  measurements.  A  transformation  of  variables  can  be  avoided  if 
the  base  transmitter  is  selected  that  has  no  cycle  slips,  otherwise  it  could  be  possible 
that  all  the  single  differenced  measurements  for  an  epoch  could  be  invalidated  by  a 
cycle  slip  on  the  base  measurement  alone.  Figure  3.5  shows  a  plot  of  all  5  transmitters 
and  their  corresponding  cycle  slips.  Other  than  transmitter  4,  all  transmitters  have 
at  least  one  cycle  slip.  Thus  number  4  is  selected  as  the  base.  With  5  measurements 
at  each  epoch,  there  will  be  4  single  differenced  measurements.  Each  will  be  the 
difference  between  the  base  and  another  transmitter. 

3.2.6  Tropospheric  Error.  Although  pseudolite  signals  do  not  have  to  deal 
with  ionosperic  error,  they  do  deal  with  much  more  error  due  to  tropospheric  effects. 
It  is  much  harder  to  measure  the  tropospheric  delay  close  to  the  surface  of  the  Earth. 
Many  models  are  only  valid  for  elevation  angles  greater  than  15  degrees.  In  the  test 
case  for  this  research,  the  greatest  elevation  did  not  get  above  6  degrees.  A  new  model 
was  used  that  can  handle  any  elevation  angle  and  relies  on  temperature,  pressure,  and 
relative  humidity  [8,9,11]. 
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Figure  3.5:  Plot  of  each  transmitter  and  the  cycle  slips  for 
each.  Other  than  number  4,  each  transmitter  experiences  at 
least  one  cycle  slip. 

The  tropospheric  delay  is  calculated  using  the  atmospheric  parameters,  pseudo- 
lite  positions,  and  estimated  receiver  location.  These  corrections  are  added  into  the 
measurement  residual  vector,  Az.  In  this  manner  the  tropospheric  error  is  modelled 
and  mitigated. 


3.3  Formation  of  the  Surface  Grid 

Developing  methods  to  use  predetermined  surface  data  to  aid  the  height  esti¬ 
mation  is  the  goal  of  this  research.  For  ground  based  applications,  the  test  area  can 
be  surveyed  using  precise  DGPS  techniques,  ensuring  centimeter-level  accuracy.  Once 
this  is  done,  a  grid  of  the  area  can  be  created  that  describes  the  surface  accurately 
enough  to  aid  the  least  squares  techniques  being  used.  While  formation  of  the  grid  is 
important  for  this  research,  the  actual  inclusion  of  this  information  into  the  solution 
is  actually  being  investigated.  Therefore  little  time  will  be  spent  on  methods  for  grid 
creation. 
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In  an  ideal  case,  the  surface  can  be  surveyed  in  a  grid  fashion  allowing  the 
surveyed  locations  to  be  in  a  grid  formation  to  begin  with.  In  the  road  test  conducted 
for  this  research,  no  such  survey  technique  was  possible.  Instead,  two  GPS  antennas 
were  placed  on  either  side  of  the  Locata  antenna.  This  GPS  data  was  logged  while 
data  was  collected  for  the  pseudolite  system.  This  not  only  allowed  a  surveyed  map 
of  points  to  be  made  of  the  road,  but  also  allowed  a  truth  trajectory  to  be  made  from 
the  average  of  the  two  GPS  solutions.  Again,  differential  GPS  methods  have  been 
applied  to  gain  centimeter  level  accuracy.  This  required  the  use  of  a  stationary  GPS 
reference  receiver  as  well. 

3-4  Determination  of  a  Height  Using  the  Surface  Grid 

Many  of  the  algorithms  developed  rely  on  determining  the  height  of  a  horizontal 
position  using  the  surface  topography.  Since  the  surface  information  is  contained  in 
a  grid  of  points,  interpolation  of  the  grid  must  be  done  to  find  the  height  at  some 
location  between  grid  points. 

Bilinear  interpolation  was  used  to  find  the  height  of  a  horizontal  position  using 
the  surface  grid.  This  method,  described  in  [18],  gives  the  surface  value  ( zest )  given  x 
and  y. 


Zest (x,  y)  =  (1  -  t)(  1  -  u)zid  +  t(  1  -  u)zi+ 1J  +  (1  -  t)uzitj+ 1  +  tuzi+  1J+1  (3.1) 


where  Xi  <  x  <  Xi+ \  and  yj  <  y  <  y]+\ ,  and 


t  = 


(x  -  Xi) 
{xi+l  -  Xi) 


(3.2) 


u  = 


( y  -  Vi) 


(3.3) 


(Vj+ 1  -  Vj) 

This  method  will  be  used  in  many  of  the  algorithms  developed  as  a  tool  to  calculate 
the  height  on  the  surface  grid. 
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3.5  Basic  Formation  of  ILS  Terms 

Creation  of  the  basic  vectors  and  matrices  for  this  research  is  described  below.  A 
normal  epoch-by-epoch  ILS  structure  will  be  introduced,  and  then  this  basic  structure 
will  be  modified  for  the  various  algorithms  that  are  implemented. 

The  measurements  that  will  be  used  primarily  in  this  research  are  single  differ¬ 
enced  carrier  phase  measurements.  Assuming  there  are  5  pseudolites  (labeled  PRNS 
1  through  5),  one  pseudolites  is  selected  as  the  base  transmitter  for  differencing.  The 
measurements  from  the  remaining  transmitters  will  each  be  subtracted  from  the  base 
transmitter’s  measurements.  With  PRN  1  designated  as  the  base  transmitter  the 
measurement  vector  will  appear  as 


z  =  [VC*  VC*  VC*  VC*]T  (3-4) 

where  the  single  difference  (V)  measurements  are  described  by  Equation  (2.7).  The 
tropospheric  error  is  removed  using  the  error  model  described  in  [11],  The  psueudolite 
clock  error  is  removed  due  to  Locata  synchronization,  while  the  multipath  and  mea¬ 
surement  noise  terms  are  combined  into  the  noise  term.  Under  these  assumptions,  an 
individual  measurement  equation  can  be  expressed  as 


=  X-^Vr^  +  VN1/  +  u  (3.5) 

where 


V 

=  noise  term 

\7N\2 

=  ambiguity  term  for  V</> 

> 

=  single  difference  range 
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and  Vr^2  can  be  represented  as 


Vr^2  =  y/ (xi  -  xA )2  +  (yi  -  yA )2  +  (-1  -  ^)2 

-a/(x2  -  xa)2  +  (1/2  -  ?m)2  +  (z2  ~  zA )2  (3.6) 

(3.7) 

[aq,  7/1,  Zi],[x2, 1/2,  Z2],  and  [xA,yA,  zA]  are  the  positions  of  pseudolite  1,  pseudolite  2, 
and  receiver  A,  respectively. 

The  state  vector,  x,  contains  the  positions  of  the  reciever  and  all  of  the  ambiguity 
states  for  the  entire  batch  process.  The  state  vector  corresponding  to  the  above 
measurement  vector  is 


X  =  [xA  VA  ZA  VN\2  VN\3  VN1/  VN\5]T  (3.8) 

The  right  hand  side  of  Equation  (3.5)  contains  the  differenced  range  from  the  esti¬ 
mated  receiver  position  to  two  transmitters.  Equation  (3.5)  can  be  expressed  in  the 
form  z  —  h  (x)  +  u,  where 


h a  (x)  =  A  1(y/ (xi  -  xA )2  +  ( Vi  ~  Va)2  +  (21  -  zA)2 

-\/(x2  ~  xA)2  +  (2/2  -  Va)2  +  (z2  ~  zA)2)  +  ViV)2  (3.9) 

Formulation  of  the  H  matrix  is  done  by  finding  the  partial  derivatives  of  (x)  wrt 
the  state  vector. 

This  is  shown  for  a  single  row  of  the  H  matrix  as  an  example. 
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(3.10) 


dhj  W 

dxA 


Xl-XA 


A  \  [(xi— XA)2+(yi— yA)2+(2i— -a)2]1^2 


X2  ~XA 


X  \  {(x2-XA)2  +  (y2-yA)2  +  {z2-ZA)2]1/2 

=  H-ei  +  c|} 


dhf  (x) 

dyA 


yi-VA 


A  l  [(*i-*A)2+(yi-yA)2+(2i-^A)2]1/2 


V2-VA 


A  l  [(^2— 31a)2 +  (2/2— 2/a)2 +  (22— Za)2]1^2 

=  j{-ei  +  el} 


(3.11) 


dhj  W 

dzA 


Zl-ZA 


A  \  [(xi-ia)2+(2/i-2/a)2+(zi-za)2]1/2 


Z2-ZA 


A  l  [(a:2— 31a)2  +  (2/2— '2/a)2  +  (z2— Za)2]1^2 

=  J  {-el  +  el} 


(3.12) 


Ml  (x) 

dVN\2 


(3.13) 


where 

eA  =  [4  el  el\  (3-14) 

is  the  unit  line-of-sight  vectors  pointing  from  the  mobile  receiver  to  pseudolite  j. 

Forming  these  partial  derivatives  into  the  complete  H  matrix  for  a  single  epoch 


is 


e1 

+  e2) 

1 

0 

0 

0 

e1 

+  e3) 

0 

1 

0 

0 

e1 

+  e4) 

0 

0 

1 

0 

e1 

+  e5) 

0 

0 

0 

1 

(3.15) 
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where  eJ  is  the  unit  line-of-sight  vector  between  the  mobile  receiver  and  pseudolite  j. 

The  measurement  covariance  matrix,  R,  is  a  square  matrix  with  diagonal  terms 
corresponding  to  the  noise  variance  in  each  measurement  and  off-diagonal  terms  cor¬ 
responding  to  the  cross-covariances  between  measurements.  Carrier-phase  measure¬ 
ments  are  the  only  measurements  being  considered  in  this  example.  Carrier-phase 
cross-covariances  are  assumed  to  be  half  of  the  carrier-phase  variance,  because  half  of 
the  measurements  are  in  common  due  to  differencing  [8]. 


where  a ^  =  .004  m2  [20]. 

From  the  example  supplied  it  is  important  to  note  the  number  of  states  and 
the  number  of  measurements.  Given  the  5  pseudolites,  there  are  4  single  differenced 
measurements  available.  With  3  position  states  and  one  ambiguity  state  for  each  mea¬ 
surement,  the  number  of  states  reaches  7.  It  is  not  possible  to  accurately  estimate  7 
states  from  only  4  measurements.  The  solution  to  this  problem  lies  in  the  commonal¬ 
ity  of  the  ambiguities.  Assuming  there  are  no  cycle  slips,  the  same  ambiguity  state  is 
associated  with  all  the  measurements  from  a  single  transmitter.  A  batch  ILS  process 
will  allow  measurements  from  multiple  times  to  be  used  in  estimating  the  position  of 
the  receiver  at  each  time  as  well  as  the  ambiguity  state  common  to  the  measurements 
across  different  times. 

3.5.1  Batch  processing.  For  a  batch  ILS  estimation  problem  the  above 
single-epoch  example  will  be  expanded.  Instead  of  using  the  measurements  from 
a  single  time  epoch,  measurements  from  multiple  epochs  will  be  used.  The  states 
associated  with  each  time  epoch  will  also  be  used  to  create  a  single  larger  state  vector. 
Consider  using  the  same  number  of  pseudolites  as  in  the  previous  section,  but  now 
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instead  of  using  only  a  single  time,  three  times  will  be  used.  The  measurement  and 
state  vectors  will  now  appear  as 


z  =  [V,y2  V 4>™  V4>a  vy ?  V*2  V*3  V*5  V0(12  V4>)i  V4>H  V*']T 

(3.17) 

x  =  [xtl  yn  zn  x„  ya  za  xa  ya  za  V.V12  V1V13  ViV14  V/V15]T  (3.18) 

where  subscripts  t±,  £2,  and  t%  indicate  the  time  epochs  for  the  particular  measurement 
or  state,  ft  is  assumed  that  there  are  no  cycle  slips,  so  the  ambiguities  remain  valid 
for  all  three  epochs.  It  is  noted  that  instead  of  estimating  7  states  with  4  measure¬ 
ments  (as  in  the  single  epoch  example),  now  13  states  are  being  estimated  with  12 
measurements.  If  more  epochs  were  considered  in  the  batch  process,  the  number  of 
states  will  be  smaller  than  the  number  of  measurements,  and  accurate  estimation  of 
the  states  can  be  accomplished. 

The  batch  H  matrix  is  a  block  diagonal  combination  of  single-epoch  H  matrices 
(from  Equation  3.15),  with  ambiguity  terms  (l’s)  consolidated  in  the  last  4  columns, 
as  shown  below 


48 


K-^i  - 

'  eu) 

[0,0,0] 

[0,0,0] 

1 

0 

0 

0 

x(-eti  - 

'  etl) 

[0,0,0] 

[0,0,0] 

0 

1 

0 

0 

x(-eti  - 

'  etl) 

[0,0,0] 

[0,0,0] 

0 

0 

1 

0 

x(-eu  - 

■e5n) 

[0,0,0] 

[0,0,0] 

0 

0 

0 

1 

[0,0,0] 

+ 

e?2) 

[0,0,0] 

1 

0 

0 

0 

[0,0,0] 

X  (  e*2 

+ 

e!) 

[0,0,0] 

0 

1 

0 

0 

[0,0,0] 

X  (  e*2 

+ 

4) 

[0,0,0] 

0 

0 

1 

0 

[0,0,0] 

iv  (  et2 

+ 

eg,) 

[0,0,0] 

0 

0 

0 

1 

[0,0,0] 

[0,0,0] 

+ 

eg,) 

1 

0 

0 

0 

[0,0,0] 

[0,0,0] 

+ 

eg,) 

0 

1 

0 

0 

[0,0,0] 

[0,0,0] 

+ 

et3) 

0 

0 

1 

0 

[0,0,0] 

[0,0,0] 

i(-4 

+ 

eg,) 

0 

0 

0 

1 

(3.19) 


The  batch  R  matrix  is  also  a  block  diagonal  combination  of  single-epoch  R 
matrices  as  shown  below 


^TdTT00000000 

f  aj  I  f  00000000 

a~i  *4  al  oooooooo 

a~i  a~i  a~i  oooooooo 

0000  <rj  ^^^0000 

0000  al  §^0000 

oooo  °i  al  ^oooo 

0000  0000 
00000000  00  ^  ^  ^ 

oooooooo  al  °i 
oooooooo  °i  ^  °i 

00000000  I  I  f  ffj 


(3.20) 
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3.6  Algorithms 

Once  the  data  has  been  analyzed  for  cycle  slips  using  the  SNR  measurements, 
the  ambiguity  states  can  be  associated  with  the  measurements.  The  initial  trajec¬ 
tory  is  calculated  using  the  pseudorange  measurements.  After  this  is  done,  the  epoch 
selection  process  is  done  to  determine  which  times  will  be  used  to  form  the  batch 
process.  Next,  the  batch  ILS  process  is  used  to  estimate  the  ambiguity  states.  During 
the  ILS  estimation,  the  tropospheric  delay  is  modelled  and  subtracted  from  the  mea¬ 
surements.  Finally,  once  the  ambiguities  are  known,  the  final  solution  can  be  found 
epoch  by  epoch  using  the  ILS  algorithm. 

Figure  3.6  illustrates  the  overall  process  described  above  in  a  block  diagram. 
This  process  is  followed  for  each  height  aiding  algorithm.  The  process  is  the  same  for 
using  either  true  data  or  simulated  data.  With  true  data,  the  surface  grid  is  computed 
from  the  GPS  truth  positions,  while  in  simulated  data  the  grid  is  set.  The  simulated 
data  is  also  generated  in  the  ENU  frame  so  that  any  frame  rotations  required  for  real 
data  are  not  needed. 

In  common  practice,  the  solution  is  constrained  to  the  horizontal  plane.  A 
nominal  height  is  assumed  and  the  problem  is  now  reduced  to  finding  the  horizontal 
positions.  In  most  cases  for  which  height  constraints  are  used,  the  height  is  known. 
In  this  research  the  receiver  is  travelling  on  a  known  surface.  The  height  is  therefore 
known  if  the  horizontal  positions  are  known.  The  following  sections  will  provide  the 
presentation  of  the  five  algorithms  developed  to  use  the  known  topography  to  aid  the 
ILS  estimation  in  this  research. 

3.6.1  Method  1:  Horizontal  Solution  Constrained  to  the  Surface.  The  first 
method  will  assume  an  initial  position,  including  an  initial  height.  This  solution 
will  be  used  as  the  linearization  point  for  the  H  matrix.  The  problem  will  then  be 
constrained  to  the  horizontal  plane  by  removing  the  height  from  the  state  vector 
and  removing  the  corresponding  column  from  the  H  matrix.  It  is  important  to  note 
that  the  full  3-d  position  will  be  used  to  calculate  the  measurement  residuals  Az. 


50 


Figure  3.6:  Block  diagram  of  overall  algorithm  process. 


In  this  manner  the  initial  height  will  not  be  estimated  in  the  ILS  process.  The 
updates  to  the  states  will  be  calculated  and  used  to  create  the  new  position  for  the 
beginning  of  the  next  iteration.  Once  the  new  states  are  determined  and  before  the 
next  iteration  starts,  the  horizontal  components  of  the  position  will  be  used  to  find 
their  corresponding  height  on  the  grid.  This  height  will  be  determined  by  using  the 
known  topography  grid. 

The  horizontal  components  are  used  to  find  the  four  surrounding  grid  points. 
Bilinear  interpolation  is  then  used  to  find  the  height  that  matches  the  horizontal 
components.  Once  this  new  height  is  found,  the  process  starts  again.  In  this  manner, 
the  solution  will  always  lie  on  the  surface. 

This  method  estimates  only  the  horizontal  components  of  the  receiver  position. 
This  is  important  in  considering  the  number  of  states  that  will  be  estimated  and 
thus  the  size  of  the  matrix  inversion  required  in  the  ILS  process.  Each  epoch  that 
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Figure  3.7:  Illustration  of  Method  1. 

is  included  in  the  batch  process  will  contribute  2  position  states.  While  the  number 
of  position  states  estimated  will  change  with  each  method,  the  number  of  ambiguity 
states  estimated  will  remain  the  same  between  the  different  methods. 

An  illustration  of  how  Method  1  might  converge  to  a  solution  is  given  in  Figure 
3.7.  The  iteration  begins  with  an  initial  position  (green  star).  The  height  is  held 
constant  at  the  initial  height,  while  the  horizontal  position  is  estimated.  Once  the 
estimation  process  has  updated  the  initial  position  (only  the  horizontal  points),  the 
solution  is  constrained  to  lie  on  the  surface  using  bilinear  interpolation.  Once  the  new 
height  is  found  using  the  surface  grid  the  process  begins  again.  This  is  done  until  a 
solution  is  converged  upon. 

3.6.2  Method  2:  3-D  solution  with  pseudo-measurement.  The  second  algo¬ 
rithm  tested  more  directly  follows  the  common  height  constraint  methods  [17,23,27]. 
A  pseudo-measurement  of  the  height  is  created  based  upon  the  bilinear  interpolation 
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of  the  current  location  on  the  grid  (Equation  (3.1)).  The  measurement  is  incorporated 
alongside  the  single-differenced  phase  measurements.  This  allows  the  solution  to  be 
found  in  3-dimensions.  A  variance  is  associated  with  the  height  measurement  and 
is  augmented  to  the  existing  measurement  noise  covariance.  This  variance  is  chosen 
to  be  (.001m)2.  This  assumes  the  grid  is  accurate  to  1  mm.  A  1  mm  error  is  ac¬ 
ceptable,  since  the  truth  height  will  also  be  fitted  to  the  surface  grid  using  the  same 
interpolation  method.  Each  position  calculated  during  the  iteration  process  is  used 
to  create  a  pseudo-measurement  of  the  height  based  on  the  pre-formed  grid.  This 
pseudo-measurement  will  be  used  in  the  next  iteration.  Instead  of  limiting  the  ILS 
estimate  to  only  the  horizontal  positions,  it  is  now  able  to  estimate  all  3  position 
states. 

Estimating  3  position  states  instead  of  2  will  increase  the  size  of  the  required 
matrix  by  roughly  50  percent.  For  instance,  a  batch  process  of  100  epochs,  with  a 
total  of  10  ambiguity  states,  will  invoke  a  210  by  210  matrix  inversion  using  method 
1,  and  a  310  by  310  matrix  inversion  using  Method  2. 

Given  the  same  conditions  described  in  previous  sections,  but  with  a  pseudo¬ 
measurement  now  available,  the  H  matrix  will  contain  an  additional  row  correspond¬ 
ing  to  the  new  measurement.  Since  the  height  measurement  is  a  direct  measure  of 
the  height  component,  its  partial  derivatives  wrt  any  other  states  are  zero.  Thus,  if 
the  state  vector  contains  the  three  positions  states,  and  four  ambiguity  states  the  row 
appended  to  the  H  matrix  will  be  [0  0  1  0  0  0  0].  A  new  diagonal  entry  will  be  added 
to  the  measurement  noise  matrix,  R,  as  well: 


R 


method2 


R 


0 


cr 


HC 


(3.21) 


where  a  value  of  <j2hc  =  (0.001m)2  was  chosen  since  the  surface  grid  was  created  using 
DGPS. 
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Figure  3.8:  Illustration  of  Method  2. 

Again,  an  illustration  will  be  used  to  show  how  Method  2  behaves  (Figure  ??). 
This  algorithm  does  not  require  that  the  solution  lies  on  the  surface.  Before  estimation 
occurs,  the  horizontal  positions  of  the  initial  point  are  used  to  measure  the  height  of 
the  surface  at  that  point  using  bilinear  interpolation.  This  pseudo-measurement  is 
incorporated  into  the  solution  process.  The  height  constraint  really  attempts  to  force 
the  solution  to  stay  near  the  height  of  the  surface  known  from  the  last  iteration.  The 
new  solution  is  estimated  in  three  dimensions. 

3.6.3  Method  3:  Single  Value  Decomposition  Epoch-by- Epoch  Transformation. 

Problems  can  occur  in  ILS  when  the  rank  of  H  is  less  than  the  number  of  states. 
This  will  cause  the  matrix  H7  H  to  be  singular  and  uninvertible. 

Another  problem  occurs  when  a  state  parameter  is  nearly  unobservable.  While 
there  is  no  rank  deficiency  causing  the  problem  of  uninvertiblity,  the  estimate  of  the 
unobservable  state  will  be  poor.  To  eliminate  this  problem,  the  state  space  of  the 
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system  can  be  transformed  to  eliminate  the  problem  states.  The  transformation  will 
change  both  the  state  vector  and  the  H  matrix  as  shown. 

Xn'xl  =  ^n'xnXnxl  (3.22) 

^mxn'  =  Hmxn  (T  ) nxn'  (3.23) 

where  n  <  n  and  the  prime  notation  indicate  the  new  state  vector  and  H  matrix. 
The  above  and  following  derivation  closely  follow  the  development  given  in  [21].  It  is 
also  very  similar  to  the  concept  of  primary  component  analysis  [12], 

3.6.3. 1  Transformation  matrix.  Determining  the  transformation  ma¬ 
trix  to  use  is  accomplished  by  exploiting  the  method  of  SVD.  First,  the  singular  value 
decomposition  of  H  is  generated 


H  =  USVJ 


(3.24) 


where 


U  = 


T  T 

Ui  U2 

I  I 


T 

ID 

I 


(3.25) 


in  which  u^s  are  left  singular  vectors  of  H  (the  right  eigenvectors  of  HH3 ), 


T  T  T 

Vi  V2  •  • ■  vn 

'LX  X 


(3.26) 
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where  v,’s  are  right  singular  vectors  of  H  (the  right  eigenvectors  of  H7  H),  and 


£  = 


O’  1 


^2 


(J  n 


(3.27) 


where  0\  >  02  >  •  •  •  >  an  >  0  are  the  singular  values  of  H  (the  square  roots  of  the 
eigenvalues  of  H7  H). 

Consider  the  H  as  a  mapping  function  that  maps  the  state  vector  into  the 
measurement  vector.  The  SVD  of  H  gives  an  orthonormal  basis  for  IT7  H  (the  right 
singular  vectors,  v,;’s  )  and  a  measure  of  the  observability  of  each  basis  vector  (the 
singular  values,  ads).  Small  singular  values  are  indicators  of  unobservability.  Creating 
the  transformation  matrix  is  as  simple  as  eliminating  the  right  singular  vectors  that 
correspond  to  the  deficient  singular  values.  As  an  example,  suppose  there  are  n 
singular  values,  selecting  k  of  the  largest  singular  values,  where  k  <  n  would  yield  a 
transformation  matrix  of 

(3.28) 

Now,  this  T  matrix  can  be  used  to  transform  the  H  matrix  and  the  state  vector  as 
shown  in  Equations  (3.22)  and  (3.23).  Using  the  new  H  matrix  and  state  vector, 
perform  the  same  1LS  method.  Once  the  transformed  states  have  been  estimated  the 
same  transformation  matrix  can  be  used  to  find  the  original  state  vector. 


x  = 


T7x' 


(3.29) 
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When  transforming  back  to  the  original  state  space,  the  estimate  of  the  states  that 
were  eliminated  due  to  near  unobservability  will  still  be  poor.  This  creates  the  need 
for  these  states  to  be  estimated  using  the  known  surface  grid,  just  as  in  Method  1. 

Method  3  will  use  the  above  SVD  approach  to  change  the  states  of  the  batch 
process  through  an  epoch-by-epoch  formulation.  As  previously  mentioned,  the  batch 
process  can  be  conceptualized  as  a  combination  of  single  epochs  into  one  ILS  process. 
The  SVD  approach  can  be  used  at  each  epoch  to  transform  the  block  diagonal  com¬ 
ponents  of  the  batch  H  matrix.  These  block  diagonal  components  are  initially  related 
to  the  three  ENU  position  states  of  a  single  epoch.  The  SVD  approach  will  transform 
the  states  into  the  principal  directions  and  eliminate  the  least  observable  state.  For 
a  ground-based  system  this  is  usually  the  vertical  state.  It  is  assumed  in  this  method 
that  there  is  a  deficiency  in  the  vertical  direction. 

At  each  epoch,  the  SVD  approach  will  be  used  to  transform  the  states  for  that 
epoch  alone.  This  will  be  done  for  each  epoch  included  in  the  batch  process.  After 
this  is  done,  the  batch  H  matrix  is  created  using  the  H  matrix  from  each  epoch. 

In  this  method  the  lowest  singular  value  will  be  selected  and  its  correspond¬ 
ing  state  will  be  eliminated.  This  will  always  constrain  the  solution  to  the  optimal 
plane  determined  by  the  SVD  approach.  Accordingly,  the  size  of  the  matrix  inversion 
required  in  the  batch  ILS  process  is  the  same  size  as  in  Method  2. 

Figure  ??  shows  how  Method  3  constrains  the  solution  to  lie  on  a  plane  (similar 
to  Method  1).  The  difference  is  that  the  plane  in  this  method  is  formed  to  eliminate 
the  least  observable  direction.  The  SVD  approach  will  transform  the  state  space  into 
a  2-D  plane  and  continue  to  iterate  similar  to  Method  1.  The  difference  only  lies  in 
the  alignment  of  the  planes. 

3.6.4  Method  4'-  Single  Value  Decomposition  Batch  Transformation.  Instead 
of  using  the  SVD  approach  on  each  individual  epoch  of  the  batch  process,  it  is  possible 
to  transform  the  entire  state  vector  and  H  matrix.  The  same  method  of  SVD  is  applied 
as  in  the  epoch  by  epoch  approach  above.  It  is  important  to  note  that  in  this  case 
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Figure  3.9:  Illustration  of  Method  3. 

a  threshold  is  set  to  determine  which  singular  values  and  thus  which  singular  vectors 
will  be  used  to  form  the  transformation  matrix.  The  threshold  is  selected  so  that  the 
solution  does  not  have  invertibility  issues.  This  means  the  third  dimension  will  not 
be  automatically  removed  unless  the  singular  value  associated  with  it  is  below  the 
threshold. 

The  surface  information  will  be  applied  at  the  end  of  each  iteration  using  the 
same  bilinear  interpolation  method  as  in  the  previous  methods.  Once  the  states 
have  been  transformed  back  to  the  original  ENU  representation  the  height  can  be 
constrained  to  lie  on  the  surface  using  the  two  horizontal  components. 

This  method  differs  from  the  previous  SVD  method  in  its  implementation  on 
the  entire  batch  process  at  once  and  its  use  of  a  threshold.  The  threshold  will  allow 
the  SVD  process  to  only  eliminate  states  that  have  singular  values  below  a  threshold 
and  more  likely  to  have  true  observability  deficiencies.  Considering  the  entire  batch 
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process  at  once  can  allow  the  observability  of  states  to  depend  upon  information  from 
other  epochs  other  than  their  own  associated  time.  Since  there  is  no  constraint  on 
the  number  of  states  that  are  eliminated,  there  is  no  set  size  to  the  matrix  inversion 
in  the  1LS  calculation,  other  than  the  maximum  number  being  3  states  per  epoch  and 
the  original  number  of  ambiguity  states. 

Since  this  method  does  not  force  the  constraining  of  the  solution  to  a  plane,  there 
cannot  be  an  easy  visualization  of  its  solution  process.  This  method  also  transforms 
all  the  position  states  from  all  the  epochs  at  once.  This  could  allow  the  constraint 
to  be  much  different  than  before.  In  essence  the  transformation  will  always  align  the 
new  state  space  with  the  directions  of  most  observability. 

3.6.5  Method  5:  Tangent  Plane  Transformation.  This  method  takes  advan¬ 
tage  of  the  topographical  surface  in  a  different  manner  than  the  previous  methods. 
Instead  of  constraining  the  height  of  the  solution  or  finding  the  optimal  plane  de¬ 
pendent  of  the  location  of  the  pseudolites  and  receiver,  this  approach  constrains  the 
solution  to  a  plane  tangent  to  the  known  surface. 

Finding  the  plane  tangent  to  the  surface  grid  is  accomplished  by  finding  a  basis 
for  the  plane  using  3  points.  The  current  solution  position  will  be  used  as  the  origin 
of  the  basis.  A  second  point  is  found  by  moving  a  small  distance  in  the  east  direction 
and  finding  the  associated  height  using  bilinear  interpolation.  Another  point  is  found 
by  making  a  small  change  in  the  north  direction  and  finding  its  corresponding  height. 
Using  the  center  point  and  each  of  the  two  calculated  points,  two  orthogonal  vectors 
are  created.  Next  the  unit  vectors  in  these  directions  are  found  and  used  as  the  basis 
vectors  for  the  new  solution  plane.  Since  these  vectors  are  orthogonal,  they  can  be 
defined  as  the  vectors  [1,0,0]  and  [0,1,0]  in  our  new  frame.  A  transformation  matrix 
into  this  new  plane  is  then  created  by  representing  these  basis  vectors  in  the  original 
coordinate  frame.  This  method  of  coordinate  frame  transformation  is  described  in  [13]. 

The  state  vector  and  the  H  matrix  are  transformed  similarly  to  Method  3  for 
each  epoch.  The  solution  is  constrained  to  lie  on  the  tangent  plane.  Once  the  solution 
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is  found,  it  must  be  transformed  back  to  the  original  ENU  frame.  As  in  the  SVD 
methods,  the  transformation  cannot  create  a  good  estimate  of  the  height  without 
additional  information.  Therefore  the  height  is  still  constrained  to  lie  on  the  surface 
grid  before  the  next  iteration  begins. 

This  method  constrains  the  solution  to  two  dimensions.  This  will  allow  the  size 
of  the  required  matrix  inversion  in  the  ILS  process  to  be  roughly  33  percent  smaller 
than  a  normal  3  dimensional  solution  method. 

The  illustration  of  this  method  is  shown  in  Figure  ??.  This  method  constrains 
the  solution  to  a  plane  as  in  Methods  1  and  3,  but  the  plane  is  now  found  to  be  tangent 
to  the  known  surface.  This  aligns  the  state  space  with  the  direction  of  slope  of  the 
terrain.  The  iteration  process  is  similar  to  Methods  1  and  3.  An  initial  position  is 
considered.  This  position  is  used  to  find  the  tangent  plane  and  transformation  matrix. 
The  transformation  is  used  to  constrain  the  state  space  to  the  plane.  The  solution 
is  found  and  then  constrained  again  to  the  surface  using  bilinear  interpolation.  Once 
this  has  occurred,  a  new  iteration  can  continue. 

3.1  Test  Setups 

The  following  section  will  detail  the  setup  for  the  real  world  test  as  well  as  the 
different  simulations. 

3.7.1  Real  Road  Test.  A  real  test  using  five  LocataLites  was  conducted 
at  the  Locata  test  facility.  The  locations  for  the  LocataLites  are  given  in  ECEF 
coordinates  in  Table  3.1.  These  coordinates  were  surveyed  using  Leica  System  500 
and  1200  GPS  receivers  and  processed  using  Leica  GEO  Office  [1].  A  2-D  map  looking 
down  at  the  horizontal  positions  of  the  LocataLites  and  the  road  is  shown  in  Figure 
3.11. 

The  LocataLites  surrounded  the  road  on  which  the  test  was  conducted.  This 
setup  had  been  created  by  the  Locata  engineers  for  the  purpose  of  conducting  tests 
on  the  road.  A  test  vehicle  was  configured  with  a  Locata  receiver  and  two  GPS 
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Figure  3.10:  Illustration  of  Method  5. 

antennas.  The  GPS  antennas  were  mounted  on  either  side  of  the  Locata  antenna. 
This  would  allow  the  truth  position  of  the  Locata  receiver  to  be  the  mean  of  the  two 
GPS  positions.  The  GPS  positions  would  be  calculated  using  DGPS  corrections  from 
a  fixed  reference  receiver. 

The  GPS  receivers  would  not  only  provide  a  truth  trajectory  for  comparison  of 
the  locata  solutions,  but  would  also  be  used  to  map  the  road  for  creation  of  the  surface 
grid.  This  required  multiple  passes  up  and  down  the  road  so  that  a  sufficient  number 


Table  3.1:  ECEF  locations  for  the  LocataLites  used 
during  the  real  road  test. 


Name 

PRN  # 

X  (m) 

Y  (m) 

Z  (m) 

South 

1 

-4431516.461 

2635737.865 

-3743255.799 

West 

2 

-4431939.325 

2636221.837 

-3742211.816 

North 

3 

-4432228.435 

2636339.206 

-3741769.434 

Beach 

4 

-4432842.397 

2635894.700 

-3741362.460 

East 

5 

-4432209.808 

2635733.788 

-3742244.268 
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Figure  3.11:  Horizontal  view  of  the  road  test  setup  in  the  ro¬ 
tated  ENU  frame. 

of  points  were  collected  for  an  accurate  representation  of  the  topography.  Since  data 
has  been  completely  post  processed,  it  was  possible  to  collect  GPS  data  for  the  grid 
and  also  data  for  the  LocataLite  system  simultaneously.  The  primary  test  in  this 
research  is  a  meandering  path  that  traverses  up  and  down  the  length  of  the  road. 

3.7.2  Single  Hill  Simulations.  A  simple  simulation  is  created  by  placing  a 
hill  at  the  origin  of  an  ENU  frame.  The  height  and  area  of  the  hill  can  be  changed 
to  alter  the  slope  of  the  hill.  The  position  of  the  transmitters  is  also  altered  so  as 
to  create  either  good  or  poor  geometries.  A  truth  trajectory  is  chosen,  which  allows 
artificial  measurements  to  be  created.  Either  perfect  measurements  are  used  or  white- 
Gaussian  noise  is  added  to  the  measurements.  The  truth  trajectory  need  not  be  a 
realistic  target  trajectory,  since  the  batch  ILS  approach  does  not  assume  any  target 
dynamics. 

This  simulation  is  used  to  test  how  the  different  methods  iterate  to  the  final 
solutions.  Methods  will  be  compared  on  their  final  solution  errors  as  well  as  the 
number  of  iterations  involved  in  finding  the  solution.  Different  simulations  can  be 
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conducted  using  the  same  surface  grid  of  a  single  hill.  Instead  of  the  truth  positions 
all  lying  at  the  peak  of  the  hill,  the  positions  can  be  randomly  scattered  over  the 
surface  of  the  hill.  Now  starting  from  the  same  initial  position  at  the  top  of  the  hill 
for  each  case,  the  simulation  can  test  how  the  location  of  a  truth  point  effects  each 
algorithm.  Individual  simulations  will  be  discussed  with  their  corresponding  results 
in  the  next  chapter. 

3. 8  Summary 

The  details  of  the  algorithms  and  calculation  process  have  been  discussed.  The 
initial  setup  involved  and  differences  in  the  proposed  solutions  has  been  introduced. 
In  the  next  chapter  the  results  from  a  real  world  test  are  discussed.  Comparison  of 
the  different  algorithms  will  also  be  presented. 
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IV.  Analysis 


4-1  Overview 

In  this  section,  results  and  analysis  of  the  solutions  will  be  presented.  The  first 
results  shown  are  those  from  the  real  world  testing  that  has  been  done.  The 
issues  that  are  associated  with  real  measurements  will  be  pointed  out  and  discussed. 
After  this,  the  analysis  will  shift  to  simulation  data.  Here,  different  methods  can  be 
compared  by  limiting  the  error  parameters  and  choosing  the  initial  and  final  positions. 

4-2  Real  Road  Test 

The  road  test  was  conducted  as  a  proof-of-concept  of  a  pseudolite-based  system. 
The  focus  of  this  test  is  to  show  that  centimeter  level  accuracy  is  attainable  and  to 
characterize  issues  in  working  with  real  data. 

4-2.1  Solution  Results.  The  road  test  conducted  consisted  of  driving  up  and 
down  the  length  of  the  road  while  swerving  from  side  to  side.  Figure  4.1  shows  the 
total  three  dimensional  error  in  the  solution  using  method  1.  Figure  4.2  displays  the  3 
error  components  using  method  1.  There  is  a  period  of  static  testing  (the  vehicle  was 
stationary  such  that  the  position  data  was  highly  accurate)  done  at  the  beginning  of 
the  run  for  roughly  half  of  a  minute.  The  graph  shows  the  error  under  nine  centimeters 
for  the  entire  run.  For  most  of  the  test  the  error  is  under  5  centimeters.  Near  five 
minutes  into  the  test,  the  error  jumps  higher.  Relating  this  graph  to  the  Dilution  of 
Precision  plot  in  Figure  4.3,  illustrates  why  there  is  a  jump  in  the  position  errors.  The 
geometry  is  least  favorable  near  the  end  of  the  run.  Despite  the  peak  in  the  error,  the 
mean  error  is  2.3  centimeters. 

The  road  in  this  test  is  nearly  flat.  While  this  means  there  are  no  real  problems  in 
the  solution  process  with  any  of  the  methods,  it  also  means  all  of  the  methods  behave 
nearly  the  same.  There  is  little  difference  in  the  results  for  any  of  the  methods. 

Methods  1,  2,  3,  and  5  all  give  nearly  the  exact  same  solutions  (probably  due  to 
the  flat  road).  Method  4  (SVD  on  the  batch  process),  gives  the  only  different  solutions, 
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as  shown  in  Figures  4.4  and  4.5.  There  is  now  a  larger  error  in  the  middle  of  the  run 
and  the  large  error  spike  in  the  other  methods  around  5  minutes  is  lower  in  the  results 
for  Method  4.  The  mean  error  is  still  slightly  above  2  centimeters.  The  difference 
in  the  solutions  could  be  due  to  the  number  of  dimensions  used  in  method  4  for 
estimating  its  positions.  Remember  this  method  uses  the  SVD  approach  on  the  entire 
batch  process  and  only  eliminates  states  that  correspond  to  singular  values  below  a 
set  threshold.  This  differes  from  the  other  methods  which  constrain  the  solution  to 
a  lower  dimensionality.  This  inherently  throws  away  some  of  the  information  in  the 
measurements.  In  the  case  seen  here,  there  may  be  times  that  the  geometry  still  has 
decent  observability  and  constraining  the  solution  to  two  dimensions  is  not  needed. 
This  idea  is  captured  in  Method  4. 

By  breaking  the  3-D  error  for  Method  1  into  components  (see  Figure  4.2),  an 
interesting  characteristic  is  seen  in  the  east  error.  Looking  at  this  east  error,  it  can  be 
seen  to  slowly  rise  and  then  fall.  The  error  seems  to  grow  as  the  receiver  travels  down 
the  road  and  decreases  on  the  way  back.  Although  the  tropospheric  error  has  been 
decreased  by  the  model  applied,  it  appears  that  some  residual  tropospheric  error  is 
still  present.  As  suggested  in  [8,9],  the  residual  tropospheric  error  will  be  modelled  as 
a  scale  factor  on  the  original  error.  This  is  done  for  the  real  test  and  the  results  for 
method  1  are  shown  in  Figures  4.6  and  4.7.  The  scale  factor  was  estimated  through 
iteration  by  hand.  The  goal  being  to  reduce  the  error  in  the  east  direction  toward  zero. 
While  the  peak  error  is  reduced  slightly,  the  mean  error  increases  to  2.7  centimeters. 
This  indicates  that  the  error  model  is  incorrect.  Research  into  modelling  the  residual 
tropospheric  error  and  estimating  potential  Locatalite  position  errors  is  found  in  [24], 

4-2.2  Finding  Cycle  Slips.  Detecting  cycle  slips  is  crucial  for  correct  deter¬ 
mination  of  ambiguity  states.  If  incorrect  ambiguity  states  are  related  to  the  measure¬ 
ments,  the  precision  of  carrier-phase  measurements  is  greatly  decreased.  As  stated 
previously,  cycle  slips  in  GPS  are  often  very  dramatic  and  usually  easily  found.  In 
the  case  of  LocataLites,  the  errors  can  occur  in  such  a  way  as  to  create  a  cycle  drift 
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Figure  4.1:  3-D  error  for  the  road  test  using  method  1. 


Figure  4.2:  Error  components  for  the  road  test  using  method 
1. 
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Time  from  start  of  run  (min) 


Figure  4.3:  Dilution  of  Precision  plot  for  the  road  test. 

as  opposed  to  an  insatantaneous  integer  cycle  slip.  Monitoring  SNR  is  the  current 
method  of  determining  cycle  slips.  A  low  SNR  value  indicates  a  potential  cycle  slip. 
When  the  SNR  drops  below  a  certain  threshold  it  is  assumed  that  a  cycle  slip  has 
occured  and  a  new  ambiguity  state  is  created.  This  method  is  appropriate  assuming 
that  the  SNR  is  a  good  measure  of  when  cycle  slips  occur. 

Table  4.1  gives  the  estimates  of  the  ambiguities  for  the  road  test.  The  PRN 
numbers  are  associated  with  the  transmitter  that  is  differenced  with  the  base  trans¬ 
mitter.  The  base  transmitter  in  this  case  is  PRN  4  and  has  no  cycle  slips  (hence,  its 
selection  as  the  base).  The  start  and  stop  epochs  indicate  the  beginning  and  end  of 
the  ambiguities  valid  region  of  time.  Note  that  the  ambiguity  estimates  only  changes 
by  hundredths  of  a  cycle  for  ambiguities  7  through  9.  There  is  most  likely  no  cycle 
slip  causing  the  invalidation  of  the  7th  ambiguity  state  and  the  creation  of  the  8th 
and  9th  ambiguity  state.  An  SNR  measurement  below  the  threshold  has  indicated 
there  could  be  a  cycle  slip  at  some  point  in  time  causing  these  incorrect  changes  of 
ambiguity  states.  This  anomaly  is  seen  with  PRNs  2  and  5  as  well. 
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Figure  4.4:  3-D  error  for  the  road  test  using  method  4. 
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Figure  4.5: 
4. 


Error  components  for  the  road  test  using  method 
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Figure  4.6:  3-D  error  for  the  road  test  using  method  1  with 

the  residual  tropospheric  error  modelled  as  a  scale  factor. 


Figure  4.7:  Component  errors  for  the  road  test  using  method 
1  with  the  residual  tropospheric  error  modelled  as  a  scale  factor. 
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Table  4.1:  Carrier-phase  Ambiguity  estimates  for 
the  road  test. 


Ambiguity 
State  # 

PRN  # 

Start  Epoch 

End  Epoch 

Ambiguity 
Estimate  (cycles) 

1 

1 

1 

443 

127265.794 

2 

2 

1 

3296 

150088.190 

3 

3 

1 

3313 

160819.296 

4 

5 

1 

690 

164495.548 

5 

5 

692 

3869 

164495.547 

6 

1 

806 

947 

127307.689 

7 

1 

1415 

1559 

127306.864 

8 

1 

1561 

1785 

127306.881 

9 

1 

1787 

2314 

127306.863 

10 

1 

2627 

2737 

127308.228 

11 

1 

2979 

3869 

127307.947 

12 

2 

3309 

3869 

150088.292 

13 

3 

3348 

3869 

160817.820 

A  further  indication  of  the  inadequacy  of  SNR  measurements  in  indicating  a 
cycle  slip  detection  is  shown  in  the  solutions.  Originally  the  position  error  using 
method  1  appeared  as  in  Figures  4.8  and  4.9  instead  of  the  position  errors  in  Figures 
4.1  and  4.2.  There  is  a  large  spike  early  in  the  test  run.  A  marker  has  been  placed  at 
the  peak  to  show  the  epoch  that  it  occurred.  From  analyzing  the  data  it  is  apparent 
that  an  ambiguity  state  has  just  been  created  at  this  time.  Referring  to  Table  4.1, 
ambiguity  6  can  be  seen  as  becoming  valid  at  epoch  806.  This  means  the  measure¬ 
ments  from  PRN  1  that  were  not  previously  available  due  to  an  assumed  cycle  slip, 
are  now  being  used.  Including  more  measurements  should  improve  the  solution,  and 
not  cause  a  spike  in  the  errors.  Searching  for  the  answer  led  to  investigating  the  true 
measurement  residuals. 

The  true  measurement  residuals  are  found  by  using  the  truth  trajectory  ob¬ 
tained  from  differential  GPS.  Given  the  truth  trajectory  the  true  measurements  can 
be  calculated.  Differencing  these  true  measurements  with  the  actual  measurements 
removes  the  effects  of  geometry  changes  from  the  actual  measurements.  The  changes 
in  the  residuals  therefore  indicate  the  presence  of  other  errors.  Ambiguity  changes 
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will  be  easily  seen  as  well  as  tropospheric  delays.  Figure  4.10  shows  the  true  measure¬ 
ment  residuals.  Matching  the  time  of  the  error  spike  in  the  position  error  plot  with 
the  same  time  on  the  measurement  residual  plot  shows  the  cause  of  the  large  error  as 
being  an  ambiguity  error. 

The  large  error  occurs  while  the  measurement  is  still  in  a  cycle  drift,  which  can 
be  seen  as  the  large  change  in  the  measurement  residuals  over  time.  This  large  drift 
corresponds  to  transmitter  1.  At  this  time  an  ambiguity  state  is  being  incorporated, 
based  on  good  SNR  measurements.  The  problem  is  that  the  cycle  drift  is  still  occurring 
and  the  ambiguity  that  is  valid  for  later  time  periods  is  not  valid  at  this  time.  This 
causes  the  large  error  that  is  seen.  Looking  at  the  component  break  down  of  the 
errors  (Figure  4.9),  the  majority  of  the  error  spike  is  in  the  north  direction.  Since 
the  ambiguity  for  transmitter  1  is  incorrect,  measurements  from  this  pseudolite  are 
going  to  cause  an  error  in  the  solution.  Refering  to  the  locations  of  the  transmitters 
(Figure  3.11),  note  that  PRN  1  is  the  transmitter  that  has  the  best  observability  into 
the  north  direction.  The  fact  that  the  error  from  the  incorrect  ambiguity  effects  the 
north  primarily  helps  corroborate  the  theory  that  the  cycle  slip  detection  method  is 
not  identifying  the  correct  slips. 

The  SNR  measurements  are  good  while  a  cycle  slip  is  occurring.  This  is  a  very 
strong  indication  the  SNR  measurements  are  not  adequate  for  tracking  cycle  slips  for 
the  system  used  in  this  test. 

4-2.3  Grid  Errors.  In  this  research  the  surface  topography  was  created  using 
surveyed  GPS  points.  The  positions  were  gathered  while  driving  the  test  vehicle  to 
collect  Locota  receiver  data.  Efforts  were  made  to  combine  enough  GPS  tracks  to 
create  an  accurate  representation  of  the  road.  Looking  at  some  of  the  data  suggests 
there  are  errors  introduced  when  using  a  surface  topography  created  from  sparse 
GPS  points.  Figures  4.11  and  4.12  show  the  GPS  trajectories  plotted  along  with  the 
road  trajectory.  The  color  of  the  road  trajectory  has  been  altered  to  show  positions 
with  medium  and  large  error.  Large  errors  are  those  that  are  above  66  percent  of 
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Figure  4.8:  3-D  error  for  the  road  test  using  method  1  in  early 

testing.  Large  spike  corresponds  to  a  incorrectly  initialized  am¬ 
biguity. 

the  maximum  vertical  error.  Medium  errors  are  between  33  and  66  percent  of  the 
maximum  vertical  error. 

Looking  at  the  three  plots  (Figures  4.11  and  4.12),  it  appears  that  the  locations 
with  large  vertical  errors  correspond  with  areas  that  have  sparse  GPS  positions  in  the 
vicinity.  The  turn  in  Figure  4.11  is  a  known  location  where  the  trajectories  pass  very 
close  to  the  end  of  the  surface  grid.  Using  a  larger  number  of  GPS  points  that  more 
completely  cover  the  test  terrain  might  allow  the  creation  of  a  more  accurate  surface 
grid.  This  could  possibly  reduce  the  variations  in  the  vertical  direction. 

4-3  Simulation  Results 

Simulations  will  be  used  to  gain  more  insight  into  the  differences  between  the 
different  surface  constraint  methods.  The  simulations  are  created  so  that  the  iterative 
process  of  each  method  is  seen.  Unless  otherwise  specified,  the  ambiguity  search 
process  is  skipped  in  order  to  focus  on  the  impact  of  the  various  surface  constraint 
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Figure  4.9:  Error  components  for  the  road  test  using  method 
1  in  early  testing.  Large  spike  corresponds  to  a  incorrectly  ini¬ 
tialized  ambiguity. 
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Figure  4.10:  True  measurement  residuals  for  the  road  test. 
True  measurement  residuals  are  the  difference  between  the  ac¬ 
tual  measurements  and  what  the  measurements  are  calculated 
to  be  using  the  truth  trajectory. 


methods  on  positioning  solutions.  The  ambiguities  are  initialized  to  the  true  values 
(known,  because  it  is  a  simulation).  This  eliminates  any  complications  in  comparing 
the  different  methods.  Isolating  a  particular  characteristic  of  the  methods  is  the  goal 
of  each  simulation.  In  many  cases,  perfect  measurements  will  also  be  used.  This 
will  allow  any  fundamental  differences  to  be  more  easily  seen.  Adding  noise  to  the 
measurements  can  mask  some  differences  in  the  methods.  Each  simulation  will  be 
discussed  and  the  results  using  each  method  will  be  discussed  at  the  same  time.  This 
will  allow  the  differences  between  the  algorithms  to  be  easier  to  compare. 


4-3.1  Flat  Simulation.  A  simple  flat  plane  is  simulated  first.  Perfect  mea¬ 
surements  will  be  used  initially  to  demonstrate  how  the  operation  of  each  method  in 
an  environment  with  surface  dynamics  is  similar.  The  surface  grid,  truth  positions, 
and  measurements  are  all  artificially  created.  Each  method  uses  the  same  information 
for  its  calculations. 
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Figure  4.11:  Vertical  errors  plotted  along  with  grid  trajecto¬ 
ries  in  a  turn.  The  larger  points  are  medium  vertical  errors  (also 
yellow),  and  large  vertial  errors  (also  Red). 
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Figure  4.12:  The  larger  points  are  medium  vertical  errors  (also  yellow),  and  large 
vertial  errors  (also  Red). 

(a)  Section  in  the  middle  of  track 

(b)  A  second  plot  of  another  section  in  the  middle  of  the  test  run 


75 


Table  4.2:  Average  number  of  iterations  assuming  a 
flat  surface. 


Method  # 

Average  number  of  iterations 

1 

3.00 

2 

3.11 

3 

2.99 

4 

3.00 

5 

2.98 

4-  3. 1.1  Perfect  Measurements.  Initially,  no  measurement  noise  or 
tropospheric  delay  has  been  added  to  the  measurements.  The  truth  trajectory  is  made 
up  of  100  points  at  the  origin.  The  simulation  selects  100  random  poins  within  a  40  by 
40  meter  square  around  the  origin.  These  points  will  be  used  as  the  initial  points  for 
the  each  time  epoch.  The  solution  errors  for  all  five  methods  are  plotted  in  Figure  4.13. 
The  errors  are  all  very  small  (on  the  order  of  10~8  meters).  The  differences  between  the 
solutions  is  really  negligible  -  all  the  methods  have  essential  found  the  correct  solution. 
The  number  of  iterations  is  also  a  key  indicator  of  performance.  As  expected,  in  this 
flat  scenario  there  is  little  difference  between  the  various  algorithms.  The  average 
number  of  iterations  is  shown  in  Table  4.2.  All  the  methods  took  approximately  3 
iterations  to  find  the  solution. 

4. 3. 1.2  Measurement  Noise.  The  previous  simulation  was  repeated, 
but  this  time  measurement  noise  was  added  to  the  carrier-phase  measurements.  White, 
Gaussian  noise  with  a  standard  deviation  of  1/12  cycles  (1  cm  error)  is  added  to  the 
measurements.  The  results  are  shown  in  Figure  4.14.  The  error  is  now  on  the  order 
of  a  centimeter,  which  is  to  be  expected  based  upon  the  noise  applied  to  the  mea¬ 
surements.  The  methods  all  still  take  an  average  of  approximately  3  iterations  to 
calculate  the  final  solution.  The  methods  all  perform  nearly  the  same  even  with  noise 
added  to  the  measurements. 

4-3.2  Hills  Simulation.  Simulation  of  a  hill  at  the  origin  is  a  basic  test  of 
the  methods.  Several  different  simulations  can  be  created  with  the  same  surface  in- 
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Figure  4.13:  3-D  error  for  the  flat  simulation  with  uo  noise. 
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Figure  4.14:  3-D  error  for  the  flat  simulation  with  noise  added. 
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Table  4.3:  Locations  of  pseudolites  in  ENU  frame 
with  the  receiver  at  the  origin. 


PRN  # 

East  (m) 

North  (m) 

Up  (m) 

1 

-100 

100 

30.5 

2 

100 

-100 

29.5 

3 

-100 

-100 

30.0 

4 

100 

100 

30.0 

formation.  It  is  interesting  to  see  how  each  method  behaves  when  trying  to  estimate 
a  position  on  the  hill  starting  from  different  initial  positions.  Alternatively,  the  meth¬ 
ods  can  be  tested  using  the  same  initial  condition  while  trying  to  estimate  different 
positions  on  the  hill.  These  different  simulations  are  aimed  at  discovering  if  different 
methods  are  better  for  different  types  of  surfaces. 

4-3.2. 1  Single  Truth  Location/P  erect  Measurements.  The  hill  that  will 
be  used  as  surface  will  be  20  meters  tall  with  a  footprint  of  roughly  20  x  20  meters. 
The  locations  of  the  pseudolites  are  given  in  Table  4.3. 

Perfect  measurements  will  be  used  to  estimate  the  position  of  the  top  of  the 
hill.  To  see  if  any  algorithms  are  more  or  less  prone  to  poor  initial  starting  positions, 
100  different  initial  positions  will  be  simulated.  The  results  are  shown  in  Figure 
4.15.  While  the  errors  are  on  the  order  of  micrometers,  it  is  interesting  to  note  that 
method  5  has  the  lowest  final  position  errors.  Method  5  also  has  the  lowest  number  of 
iterations.  The  number  of  iterations  at  each  time  and  for  each  method  are  shown  in 
Figure  4.16.  Method  5  consistently  takes  five  iterations  to  converge,  while  the  other 
methods  take  close  to  6  and  7  iterations  on  average. 

Method  5  seems  to  be  superior  to  the  other  methods  of  aiding.  This  could  be 
due  to  its  ability  to  incorporate  the  information  contained  in  the  surface  topography 
more  efficiently.  All  of  the  other  methods  described  in  this  research  incorporate  the 
height  information  into  the  solution  process  in  a  delayed  manner. 

For  example  consider  Method  1.  The  height  is  constrained  to  lie  on  the  surface, 
but  the  height  is  found  using  the  horizontal  components  from  the  previous  iteration. 
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Figure  4.15:  3-D  error  for  the  hilltop  simulation  with  no  noise. 
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Figure  4.16:  Number  of  iterations  until  convergence  for  the 
hilltop  simulation. 
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This  does  not  effectively  constrain  the  height  of  the  current  solution  to  where  the 
solution  is  going  to  be.  As  a  result,  the  iteration  must  converge  in  steps.  In  cases  of 
near  horizontal  surfaces,  the  effect  of  this  constraint  issue  is  not  seen.  When  a  surface 
with  a  higher  slope  is  introduced,  problems  occur.  On  a  surface  with  a  high  slope,  a 
small  horizontal  change  causes  a  large  change  in  the  height.  This  can  require  more 
iterations  for  convergence  and  possibly  lead  to  a  larger  final  error. 

Method  5  is  the  only  method  that  constrains  the  solution  to  a  plane  that  is 
tangent  to  the  surface  topography.  In  the  case  of  a  surface  with  a  large  slope,  the 
solution  is  now  being  found  in  a  plane  that  is  rotated  to  align  with  the  only  possible 
direction  of  change  -  along  the  surface. 

While  the  tangent  plane  may  be  good  for  steep  slopes,  it  is  pustulated  that  it 
may  have  convergence  problems  at  the  peaks  of  hills.  If  the  peak  of  the  hill  is  very 
sharp,  a  small  change  in  the  position  would  create  a  vastly  different  solution  plane, 
causing  the  solution  to  oscillate  about  the  peak  without  converging. 

4 -3. 2. 2  Single  Initial  Position.  This  simulation  will  better  investigate 
the  different  algorithms  ability  to  estimate  positions  at  different  locations  on  the  hill. 
A  single  initial  position  is  used  for  each  time  epoch,  with  100  different  truth  positions 
scattered  over  the  hill.  To  create  a  surface  with  higher  slope,  the  hill  height  has  been 
increased  to  40  meters.  Measurements  are  still  considered  with  no  errors. 

Figure  4.17  illustrates  the  errors  in  the  final  positions  estimated  by  each  method. 
It  is  obvious  that  method  five  has  outperformed  the  others.  Looking  at  the  average 
number  of  iterations  in  Table  4.4,  there  is  not  as  strong  a  separation  in  the  number 
of  iterations  as  shown  when  estimating  the  top  of  the  hill. 

4-4  Summary 

This  section  has  shown  the  results  of  the  real  road  test  and  several  simulations. 
The  road  test  has  proven  that  centimeter  level  accurate  solutions  can  be  obtained 
using  the  algorithms  developed  in  this  research.  Analysis  of  the  real  data  has  shown 
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Table  4.4:  Average  number  of  iterations  in  hill  sim¬ 
ulation  assuming  common  initial  positions. 


Method  # 

Average  number  of  iterations 

1 

3.52 

2 

5.07 

3 

3.17 

4 

3.17 

5 

3.44 
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Figure  4.17:  3-D  error  for  the  hill  simulation  with  no  noise 

and  common  initial  positions. 
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the  need  for  better  cycle  slip  detection  methods.  Multiple  simulations  have  shown 
the  superiority  of  the  tangent  plane  constraint  method. 
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V.  Conclusion 


5. 1  Overview 

This  research  aims  at  developing  different  methods  of  aiding  and  addressing  the 
vertical  deficiencies  that  occur  in  position  estimation  using  pseudolite-based 
reference  systems. 

In  a  ground-based  reference  system  using  pseudolites,  the  vertical  deficiency  can 
be  much  greater  than  in  GPS  cases.  The  use  of  height  constraints  is  a  viable  method 
of  aiding  the  navigation  solution. 

Height  constraints  for  aiding  navigation  solutions  have  been  used  in  marine  en¬ 
vironments  in  past  research  [17,23,27].  This  is  valid  since  the  height  of  the  sea  and 
thus  the  receiver  can  be  well  known.  In  cases  of  navigating  on  a  changing  ground  to¬ 
pography,  knowledge  of  the  surface  height  can  be  attained  from  available  surface  grids 
(DTED),  or  from  surveying  the  terrain  manually.  This  knowledge  is  not  applicable 
as  in  the  marine  case.  Finding  the  height  from  a  surface  grid  requires  the  knowledge 
of  the  horizontal  position  components. 

Five  different  methods  have  been  developed  and  discussed  in  this  research. 

•  Method  1:  Constrain  the  solution  to  the  horizontal  plane  and  use  the  horizontal 
solutions  to  find  the  height  of  the  known  surface  at  those  coordinates. 

•  Method  2:  Estimate  in  three  dimensions  using  a  pseudo- measurement  to  con¬ 
straining  the  height.  The  height  constraint  will  change  during  the  iteration 
process  as  the  pseudo-measurement  is  found  by  projecting  the  nominal  solution 
onto  the  known  topography. 

•  Method  3:  Use  the  method  of  Singular  Value  Decomposition  to  remove  the  least 
observable  position  state  from  the  state  space.  This  is  done  at  every  epoch. 

•  Method  4:  Again  use  SVD  approach,  but  instead  use  it  on  the  entire  batch 
process.  This  will  remove  any  states  that  are  poorly  observable. 
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•  Method  5:  The  solution  is  constrained  to  lie  on  a  plane  that  is  tangent  to  the 
known  surface  topography.  This  plane  will  be  tangent  at  the  current  solution 
and  will  update  during  the  iteration  process. 

These  methods  were  used  to  estimate  the  position  of  a  receiver  in  a  real  road 
test  as  well  as  in  simulation. 

5.2  Conclusions 

The  road  test  was  conducted  as  a  proof  of  concept  of  using  a  pseudolite-based 
sysem  as  well  as  for  generation  of  real  measurement  data.  This  data  can  be  used 
for  the  estimation  of  real  error  parameters  such  as  residual  tropospheric  error  and 
transmitter  position  errors.  The  research  used  a  batch  ILS  estimation  process  to  cal¬ 
culate  the  floating  point  ambiguities.  Using  the  above  methods  it  has  been  shown 
that  sub-decimeter-level  accuracy  is  attainable  with  the  current  system.  Future  devel¬ 
opments  in  both  Locata  hardware  and  integer  ambiguity  resolution  algorithms  have 
the  potential  for  even  higher  levels  of  accuracy. 

The  cycle  slip  detection  method  of  using  SNR  measurements  has  been  analyzed 
and  found  lacking.  The  shortcomings  of  the  this  method  have  been  shown  decisively. 
Tropospheric  correction  has  been  done  using  a  new  algorithm  that  can  handle  low 
elevation  transmitters.  The  residual  tropospheric  delay  is  present,  but  has  not  been 
a  focus  of  modelling  in  this  research. 

The  methods  for  aiding  the  solution  process  using  knowledge  of  the  surface 
topography  were  developed  and  then  compared  using  multiple  simulations.  A  flat 
horizontal  surface  allows  all  the  methods  to  behave  the  same.  Their  identical  perfor¬ 
mance  in  the  presence  of  perfect  measurements  and  noisy  measurements  illustrates 
the  idea  that  these  methods  are  not  needed  for  a  case  when  the  height  of  the  receiver 
is  known  accurately  (as  in  marine  cases). 
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Various  simulations  using  a  single  hill  were  conducted.  The  effect  of  a  higher 
sloped  surface  and  the  addition  of  noisy  measurements  indicated  the  inadequacy  of 
using  a  pseudo- measurement  for  incorporating  the  height  constraint. 

Using  a  simulation  focused  on  attaining  solutions  at  different  locations  on  a  hill, 
it  was  found  that  constraining  the  solution  to  a  plane  tangent  to  the  surface  of  the 
hill  allowed  a  more  accurate  solution  to  be  found. 

While  the  simulations  conducted  barely  scratch  the  surface  of  possible  tests  that 
can  be  conducted,  the  method  of  constraining  the  solution  to  lie  on  a  tangent  plane  to 
the  known  surface  has  shown  the  best  convergence  in  terms  of  accuracy  and  number 
of  iterations.  This  method  incorporates  the  knowledge  of  the  surface  in  the  most 
efficient  manner. 

5. 3  Contributions 

The  following  contributions  have  been  made  by  this  research. 

•  The  largest  contribution  of  this  research  is  the  development  of  the  5  differ¬ 
ent  methods  for  aiding  the  height  determination  in  pseudolite-based  reference 
systems.  These  methods  can  be  used  to  take  advantage  of  the  known  surface 
topography  of  a  test  area. 

•  Development  of  a  batch  ILS  approach  to  determine  absolute  position  using  a 
pseudolite  network  and  carrier-phase  measurements. 

•  Characterizing  the  errors  involved  in  cycle  slip  detection  of  LocataLite  mea¬ 
surements  from  SNR  values.  Using  post  processed  data  with  a  known  truth 
trajectory  it  can  be  shown  that  this  method  is  not  sufficient  for  detection  of 
cycle  slips,  which  is  crucial  for  attaining  precise  solutions. 

•  Development  of  various  simulations  for  comparing  different  methods  of  height 
aiding.  These  simulations  are  important  for  understanding  the  tradeoffs  in  using 
each  method. 
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•  Determination  of  the  tangent  plane  transformation  method  as  the  best  method 
for  aiding  the  height  determination  in  a  pseudolite-based  reference  system. 

5-4  Recommendations 

This  research  focused  on  aiding  a  batch  least-squares  estimation  algorithm  in 
determining  the  position  of  a  receiver  using  a  pseudolite-based  reference  system.  The 
need  for  such  a  pseudolite-based  reference  system  is  great  and  further  research  is 
needed  to  bring  this  into  final  stages  of  development.  Recommendations  for  further 
research  are  given  below. 

•  Further  investigation  into  the  strengths  and  weaknesses  of  the  different  methods 
developed  in  this  research.  This  motivates  the  optimization  of  the  algorithm 
code  so  accurate  comparison  of  convergence  times  can  be  done. 

•  Implement  the  methods  developed  in  this  research  in  real  world  tests  that  con¬ 
tain  more  varied  surfaces.  The  road  used  in  this  test  was  nearly  flat  and  therefore 
all  the  methods  performed  similarly.  A  more  dynamic  test  surface  will  better 
explore  the  differences  in  the  methods  in  a  real  world  environment. 

•  More  investigation  of  the  ability  of  the  batch  least-squares  process  in  estimating 
system  errors  such  as  residual  tropospheric  error  and  pseudolite  position  errors. 
The  focus  of  this  research  was  on  development  of  the  different  surface  constraint 
methods.  More  potential  lies  in  using  the  batch  process  to  estimate  some  of  the 
system  errors  that  might  be  more  observable  using  a  batch  process. 

•  Instead  of  creating  a  surface  grid  using  GPS  survey  points,  investigate  other 
methods  of  generating  a  map  of  the  topography.  Using  Digital  Terrain  Elevation 
Data  (DTED)  information  is  an  obvious  example  that  could  be  researched. 

•  Investigate  the  use  of  implementing  these  method  for  use  in  cases  other  than 
pseudolite-based  systems.  In  marine  situations,  height  constraints  are  com¬ 
monly  used,  since  the  height  is  well  known.  The  use  of  aiding  a  GPS  solution 


with  a  known  topography  surface  (DTED  for  example)  could  be  used  to  boost 
accuracy  and  robustness  issues  when  GPS  geometry  is  poor. 

•  Using  multiple  methods  in  the  search  process  may  allow  for  more  accurate  solu¬ 
tions.  Given  a  trajectory  that  contains  areas  of  high  and  low  surface  dynamics, 
different  methods  may  be  converge  quicker  or  attain  better  accuracy  in  differ¬ 
ent  locations.  Combining  methods  could  take  advantage  of  the  strengths  and 
weaknesses  of  each  aiding  method. 
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